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SUMMARY 


The NLPAN computer code uses a finite-strip approach to the analysis of thin-walled prismatic 
composite structures such as stiffened panels. The code can model in-plane axial loading, transverse 
pressure loading, and constant through-the-thickness thermal loading, and can account for shape 
imperfections. The NLPAN code represents an attempt to extend the buckling analysis of the VIPASA 
computer code into the geometrically nonlinear regime. Buckling mode shapes generated using VIPASA 
are used in NLPAN as global functions for representing displacements in the nonlinear regime. While 
the NLPAN analysis is approximate in nature, it is computationally economical in comparison with 
finite-element analysis, and is thus suitable for use in preliminary design and design optimization. This 
document provides a comprehensive description of the theoretical approach of NLPAN, highlighting the 
capabilities developed under NASA Grant NAG-1-1215. A discussion of some operational consider- 
ations for the NLPAN code is included. NLPAN is applied to several test problems in order to demon- 
strate new program capabilities, and to assess the accuracy of the code in modelling various types of 
loading and response. User instructions for the NLPAN computer program are provided, including a 
detailed description of the input requirements, and example input files for two stiffened-panel config- 
urations. 
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1. INTRODUCTION 


Geometrically nonlinear analysis is generally required to accurately predict the ultimate strength 
of compressively loaded thin-walled structures such as stiffened panels. The nonlinear analysis of such 
configurations using finite elements is expensive in terms of computer resources. This has spurred the 
development of finite strip methods [1-11] which are less general than the Finite element method, but 
which are computationally economical, and thus are suitable for use in iterative applications such as 
preliminary design and design optimization. These methods vary widely with regard to the specific plate 
or shell theory used, the geometric modelling flexibility, whether or not imperfections are accounted for, 
the types of shape functions used, the targeted modes of response, the mathematical formulation used, 
and the solution strategies employed. No attempt is made here to review or compare all of these 
methods. A discussion of some of these methods can be found in [7], 

This document concerns one such method which is implemented in the computer code NLPAN. 
In philosophy, NLPAN represents an attempt to extend the buckling analysis of the VIPASA code [12] 
(the analysis code within the PASCO stiffened-panel design code [13]) into the postbuckling and ge- 
ometrically nonlinear regime. NLPAN uses buckling mode shapes obtained from VIPASA as global 
shape functions for representing displacements in a geometrically nonlinear analysis. NLPAN uses a 
stationary total potential energy formulation to obtain a set of nonlinear algebraic equations governing 
equilibrium. These equations have load-independent coefficients and a relatively small number of vari- 
able modal amplitudes, allowing rapid exploration of the nonlinear regime. For simple rectangular plate 
configurations, NLPAN degenerates to an exact series solution method for the von Kaiman nonlinear 
plate equations. 

NLPAN is limited to the treatment of prismatic structures which can be represented as assemblages 
of rigidly linked plate strips. NLPAN can model the static elastic response of a structure to in-plane 
normal loads, transverse pressure loads, and thermal loading which is uniform through the thickness of 
each plate strip. Figure 1 depicts a typical configuration with the applied loadings. Shape imperfections 
can be simulated. The side edges and longitudinal ends may have any of a variety of different support 
conditions. The individual plate strips may have orthotropic material properties, suitable for character- 
izing the elastic properties of some classes of laminated composite plates. 

The first comprehensive description of NLPAN was given in [7]. Under NASA Grant 
NAG-1-1215, a number of improvements and additions to NLPAN have been developed relative to the 
method as described in [7]. The major developments are listed here, along with the sections in this 
document (if applicable) in which each is discussed: 

1. Transverse pressure loading capability was added to the NLPAN code. 

2. Thermal loading capability was added to the NLPAN code. Temperature is assumed to be constant 
through the thickness of each plate strip, but may be different from plate strip to plate strip (Sect. 
3.2, 3.5). 

3. Advanced nonlinear solution strategies were incorporated into the NLPAN code which enable 
asynchronous application of the various load types, navigation of limit points, and navigation of 
simple or compound solution branch points (Appendix D). 

4. Methods for modelling a variety of support conditions at the longitudinal ends of a structure have 
been added to NLPAN. The new types of support include rotationally elastic support, eccentric load 
application, and clamped end support (Sect. 3.3-3.5). 

5. The NLPAN code was modified to decrease both the computer memory usage and the execution 
times for a given problem. Memory usage was decreased through the implementation of a data 
storage method which sets the dimensions of large data arrays based on the actual needs of a given 
problem (Sect 5.2). Execution times were reduced by rewriting key subroutines to run more effi- 
ciently. 
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6. An automated procedure was developed which unites NLPAN and PASCO (in which VIPASA is 
embedded), so that buckling-mode generation requests are passed automatically from NLPAN to 
PASCO, and material property data and buckling mode information are passed automatically from 
PASCO to NLPAN. (Assistance in performing this task is noted in the last paragraph of this sec- 
tion.) 

7. Troublesome aspects of the second-order displacement fields have been studied in detail, and sug- 
gestions for improving the method of computing the fields have been provided. Modifications were 
made to the method used by NLPAN in computing these fields (Sect. 4). 

8. The theory was developed for performing a transient dynamic analysis in situations where unstable 
critical points are encountered. The analysis has not been fully implemented in the NLPAN code, 
in part due to time constraints, and in part because the analytical capability is deemed to have little 
practical value within the limits of the NLPAN code (Sect. 3.1). 

9. Improvements and corrections were made to the technique used to control the relative values of the 
two generalized in-plane load components (Sect. 2.6, 2.7). 

10. A stationary total potential energy formulation was adopted in favor of the original virtual work 
formulation [7]. This was necessary in order to give the governing equations certain properties 
necessary to permit use of the advanced nonlinear solution strategies (Sect. 2.8). 

This document provides a comprehensive description of the capabilities of, and theory behind, the 
NLPAN computer code. Section 2 contains the details of the theory for the baseline method. Section 3 
documents the theory behind some of the recent additions to the method. Section 4 contains a discussion 
of various issues related to the computation of the second-order displacement fields, and describes certain 
modifications to the original method of computing them which have been implemented. Section 5 pre- 
sents a brief discussion of some operational considerations for the NLPAN code. Section 6 presents re- 
sults from applications of NLPAN which demonstrate program capabilities that are not demonstrated in 
other documents. Appendices A-C provide background material for the main text. Appendix D contains 
a description of the implementation of advanced nonlinear solution strategies within NLPAN. Appendix 
E contains detained user instructions for the NLPAN computer code. Other work related to NLPAN and 
performed under NASA Grant NAG-1-1215 is reported in [16], [23-24], 

In all work reported here, access to the buckling analysis of the VIPASA code [12] is obtained 
through the PASCO code [13], and therefore this document contains references to both codes. Which 
code is referenced is a particular discussion depends on the context; it should be kept in mind that 
PASCO is used only for the purpose of obtaining access to the embedded VIPASA analysis. 

Much of the work of linking PASCO with NLPAN was performed as a part of a separate research 
project by Ms. Christine Perry, graduate assistant to Professor Zafer Gurdal of Virginia Polytechnic In- 
stitute and State University. The NASA technical monitor for this project was Dr. James H. Starnes Jr. 
of the Aircraft Structures Branch of NASA Langley Research Center. 
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2. THEORETICAL APPROACH OF NLPAN 


2.1 Nonlinear Plate Theory 

Consider a rectangular plate strip and an associated local ( xyz ) coordinate system, where the x-y 
plane lies at the mid-surface of the undisplaced plate. The plate has dimensions L and b in the x- and 
y-directions, respectively, as shown in Figure 2. Displacements of the mid-surface are denoted by 
{«} = L«(*.y) v(*,y) w(x,y)J where the three components correspond to the x-, y-, and z-directions, re- 
spectively. The applied loads acting on the plate strip include force and moment resultants applied at 
the plate edges and uniform pressure Q acting in the local z-direction. 

2.1.1 Strain-displacement relations. The in-plane strain components {e} = [e, e, are restricted 
to be small compared to unity, but rotation angles associated with out-of-plane deflections and in-plane 
rotations are permitted to take on moderately large amplitudes. The Kirchhoff-Love assumptions are 
used to determine the distribution of the in-plane strain components through the thickness of the plate 
The in-plane strains are then given by 

{£} = {e c ) +z{k c } (211) 

where the mid-surface strains, {e c } = []e£ ej , and the mid-surface curvatures, fie'} = Tici ic K$,] r are 
given by 



The derivation of the strain-displacement relations is given in [7]. 

An imperfect plate is assumed to be free of internal stress resultants, and have non-zero mid-surface 
displacements {«} = [u°) which describe the imperfection shape. Let equation (2.1.2) define the func- 
tions {£'} = { £'({«})} and (k'J = {tc(w)}. Mid-surface "mechanical" strains {£-) and cultures ftc] 
are defined by 

{£ } = {£ c ({«})) - {e c ({« 0 })} , (k™) = (k c (w)} - {^(w")} (2.1.3) 

These are the strain and curvature measures to which stress resultants are proportional. 

2.1.2 Equilibrium equations. Force resultants (forces per unit length) /„ /„ and f are assumed to 
act along the edges of the plate in the x-, y-, and z-directions, respectively, and distributed bending 
moments M, and M, are assumed to act along the ^-normal and y-normal edges, respectively. Force 
resultant f is assumed to account for the Kirchhoff equivalent shear force due to a distributed edge 
twisting moment M^. Special notation is used to represent the load and displacement measures along 
the side (y-normal) edges. Each side edge of the strip is given an index number e (e = 1, 2), such that 
the edge station y, and the associated y-normal unit-vector component n, are given by 





[o, -1] 

lb, 1] 


for e = 1 
for e = 2 


(2.1.4) 


The generalized displacements and generalized force resultants along edge number e are denoted as 
{«.} = L«, v, w, y,J and \f,} = , respectively. The conventions for these parameters are 

shown in Figure 2. The edge rotation angle vy, is the derivative of w with respect to y. 

In addition to the edge loads, a pressure load Q, acting in the z-direction, may be present . The 
plate material is assumed to be linearly elastic. The transverse shear strains, and y„, are neglected, 

3 



consistent with the Kirchhoff-Love assumptions, and the transverse normal stress, a„ is neglected. The 
total potential energy jc is given by 


- f Q w dA - f [£u +f y v +f z w - njl x w,J \ I=0L dy 

J A J 0 

~ f [/«“* +fye v e +fu w e + "WJ \ e=:2 dX 

Jo * 


(2.1.5) 


where A is the planform area of the plate strip, hats denote specified edge-loads, and stress resultants 
{/V} = [_N,N, A/qJ and {A/} = [Af, M, M „] r are defined by 


{ } -i { s ) »■ 


z) dz 


( 2 . 1 . 6 ) 


Equations governing equilibrium are obtained by requiring stationary total potential energy: 

87c = 0 (2.1.7) 

By evaluating equation (2.1.7) and applying Green’s Theorem, the field equations governing plate equi- 


librium are found to be [7] 


N&X + y,y + (NyU^y = 0 

(2.1.8a) 

^ xy*x ^y'y + (^Yx^u) *x "" ® 

(2.1.86) 

^x'xx ^^xy*xy ^y*yy (YxH** x ^xy^*y)*x xy^ ^ ^ yW *y) *y l2 “ 0 

(2.1.8c) 

On the Jt-normal edges, where = ± 1, the expressions for the edge force resultants 
stress resultants are found to have the form 

in terms of plate 

fy ~ ^hff^xy + 

fz “ ^x(^Kr*x ^^xy'y ^ + ^xy^’y) 

(2.1.9) 


Along y-normal edges, the generalized force resultants are given by 
fz, ~ ^y(Axy 4" ^y^'y) 
fy . = "X 

4 = + M ry + + N y w ’y) 1,2 (2 ' 1 ' 1 0) 

m e = - rtyMy 

where the value n, corresponding to e is given in equation (2.1.4). 

2.1.3 Plate constitutive equations. The plate elastic properties are assumed to be those of a bal- 
anced, symmetric laminated composite plate, with the additional limitation that the bending/twisting 
coupling terms are neglected. The bending/twisting coupling terms become small in significance if each 
balanced pair of angle plies in a laminate span a portion of the laminate thickness which is small com- 
pared to the thickness of the laminate. The following plate constitutive equations then relate the plate 
stress resultants to the mid-surface mechanical strains and curvatures: 

[N] = [v4]{e m } , {A/} = [D]{K ,n } (2.1.11) 


where 
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A \\ A n 0 D u d \2 0 

[-4} = A 12 A 22 0 , [£>] = D n D 2 2 0 (2.1.12) 

0 0 A^j 0 0 Dffi 

12 Characteristics of Multiple-Plate Configurations 

A linked-plate structure has an associated global coordinate system, as shown for a sample config- 
uration in Figure 3. The global and local coordinate systems share a common x-axis. The overall width 
in the global y-direction is B. A set of "node lines" are defined which are analogous to nodes in a finite 
element model. The generalized displacements {«,} along the side edges of a plate strip are rigidly linked 
to the generalized displacements {(/") = \_U* V* W" of one of the node lines. The latter displace- 
ments are defined with respect to the global coordinate directions. The configuration shown in Figure 
3 has four node lines, which are labeled in the figure with circled numbers. 

The most general orientation of the side-edge of a plate strip relative to the associated node line is 
shown in Figure 4(a). There is an angle p between the local and global y-z coordinate directions, and 
offset measures e y and e,. Using offsets allows an accurate representation of complex cross sections, as 
can be seen in the example shown in Figure 4(b). The generalized displacements {u,} can be expressed 
in terms of {(/*} by applying sequential transformations accounting for rotation, and eccentricity, re- 
spectively (Wittrick and Williams [12]): 

u e = U n 

v e =V n cos p - W n sin p 

w e = V n sinp+W' 1 cos p (2.2.1) 

u e Ug * 4 - g z Wg, x + Cy v e , x 

v e = ig + e 2 y e 

W e = W e ~ CyVt ( 2 . 2 . 2 ) 

The second transformation above must evaluated after the functional form of the displacements has been 
specified. 

The generalized force resultants {f } acting along the side edge of a plate strip can be transformed 
to statically equivalent generalized force resultants {/"} = ] which act at the corresponding 

node line and which are defined with respect to the global coordinate directions. The transformation 
which relates [f t } and {/"} is obtained by equating the virtual work due to \f t ) acting through the virtual 
displacements {6 w,} with the virtual due to (/"} acting through the virtual displacements {5t/"}. The 
vector of generalized force resultants along a node line {Z^"} = F* F* Af"] equals the sum of the 
contributions {/*} from all plate edges terminating at the node line. Externally applied loads acting along 
a boundary node line are not included in (F"). 

Along the boundaries of the panel, it is assumed that the only non-zero specified loads are the in- 
plane normal components. The total potential energy for the entire structure can then be expressed as 
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(2.2.3) 



where p is the index number of a plate strip, P is the total number of plate strips in the structure, and 
n, and n 2 are the two designated boundary node-lines. The practical details of specifying the imposed 
loads are described in Section 2.8.2. ^ 


23 Loading and Boundary Conditions 

2.3.1 Global loads. N* is the total x-normal end load divided by the width B. N# is the total 
y-normal edge load divided by the length L. These global load measures are defined to be positive for 
tension. A non-zero load N# may only be applied across a continuous, flat skin, because the initial re- 
sponse to in-plane loading is assumed to be free of plate bending. The effective changes in the panel 
length and width associated with applied loads are denoted A u and Av, respectively, where Av is meas- 
ured between two designated boundary node lines. (In cases where the boundary node-lines deform Av 
is a mean value.) 

In a nonlinear analysis, the generalized in-plane load parameter X controls either Au or N& de- 
pending on whether input parameter CONTRL is specified as 'D* or ’L\ respectively (displacement 
control or load control). For configurations which admit biaxial loading, two options are available for 
controlling the second load component. Table 1 summarizes these options. In the table, A 14 , Av t , Njgl, 
and N^, form a set of self-consistent parameter values corresponding to a unit solution for linear! 
Un ^ U M^ e< ^ response. A unit load system is specified in terms of two parameters, one from the pair A Uc 
ana N ml, and one from the pair Av L and N^_. The two unspecified parameter values can be computed 
as described in Section 2.5. The option selected for control of the generalized in-plane loading (Table 
1 ) determines which two parameters should be specified. 

Illustrations of how the control options of Table 1 are used to simulate a variety of different sce- 
narios for generalized in-plane loading are provided in Table 2. The four different cases illustrated in 
Table 2 are: A) loading at a constant ratio NjcJN^, B) control of load component while holding Av 
to zero, C) loading at a constant ratio Av/Au, and D) uniaxial loading for configurations which do not 
admit biaxial loading. For Case D, N# is zero in prebuckling response, but a non-zero load N# may arise 
m the nonlinear regime if the boundary node-lines are restrained. 

Transverse pressure loading is permitted. The loading Q is assumed to be uniform over any given 
plate strip, but may be applied selectively to plate strips as necessary to simulate the desired loading. 
A load parameter p controls the amplitude of the pressure loading independently of the in-plane loading 
as discussed in Section 2.8.3. 

2.3.2 End-support conditions. The effective boundary support condition at the panel ends is deter- 
mined by the harmonic form of the buckling mode shapes (obtained from the buckling analysis of 
VIP ASA) used as global shape functions. The transverse displacements v and w have a sinusoidal de- 
pendence on x (see Section 2.6), and thus the buckling modes are actually those of an infinite-length 
structure which is supported against transverse displacements at uniform intervals. Axial displacement 
u has a cosine dependence on x. When the buckling results are applied to a finite-length structure, the 
ends of some panels may both rotate and warp during buckling (and postbuckling). This makes the 
concept of a neutral bending axis somewhat ill-defined; nonetheless a generalized neutral axis can be 
considered to exist, as determined by the points of zero-axial-displacement (for a buckling mode) at the 
longitudinal ends of the panel. The interpretation offered here is that each longitudinal end is loaded by 
a generalized knife-edge support which acts along the generalized neutral axis. The effective axial dis- 
placement of the panel end is then the axial displacement of the knife-edge. This issue is discussed 
further in Section 2.8.2. 
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The transverse displacements v and w of the buckling eigenfunctions are zero at the ends of the 
panel; in addition, the local bending moment M, is uniformly zero at the panel ends, implying a simple 
support condition for each plate strip. In general, however, the local edge shear-force resultant N n is 
non-zero, and can not be independently controlled. 

2.3.3 Boundary conditions along the node-lines. Two of the node lines are designated as boundary 
node-lines. The index numbers of these two node lines are denoted n, and n 2 corresponding to the panel 
edges which have edge-normal unit vector components rt, = - 1 and *,= 1, respectively, with reference 
to the global coordinate axes. Along each of the boundary node-lines, four boundary conditions are 
specified corresponding to the four degrees of freedom along a node line. In general terms, either a 
displacement condition (BCVEC=1) or a load condition (BCVEC=2) can be specified for each of the 
degrees of freedom. For the degrees of freedom corresponding to {/", W", and 4 / ", only homogeneous 
boundary conditions can be specified (zero generalized displacement or zero generalized load). For the 
degree of freedom corresponding to V*. non-homogeneous boundary conditions may be imposed con- 
sistent with the limitations given in Table 1. 

The choices for boundary conditions are expressed mathematically in Table 3. The term V L re- 
presents the node-line displacement associated with the linear response of an imperfection-free structure 
to the unit in-plane loads. A third boundary-condition option (BCVEC=3) is available for the degree 
of freedom corresponding to V" (Component 2). The symbol F; appearing in the table is the mean value 
of the force resultant F; over the length L. With this third boundary condition option, the net edge- 
normal load is controlled while the edge is constrained to remain straight with respect to in-plane dis- 
placements. The option is useful in modelling symmetry conditions, or in modelling an edge which is 
reinforced so as to remain straight. 


For the non-boundary node-lines to be in equilibrium, the generalized force-resultants must vanish: 


m = {0} 


* = 1,2 N 

n*n x ,ri2 


(2.3.1) 


2.4 Expansion of the Displacement Functions 

Displacements are assumed to have the following general form on each plate strip: 
{u)=k[u L ]+q i {u i }+q i qj{u i j} (/,;'= 1, 2, ... ) (2.4.1) 

where summation over repeated indices is implied. Symbol X is the generalized in-plane load parameter, 
and {14.} = £uc v t 0] is the displacement solution for linear, unbuckled response to specified unit, 
global, in-plane loads. The i* buckling mode shape is {«,} = £«,- v.-w,-] , and the associated "modal am- 
plitude" is q>. Shape functions {%} = [% v,, w,J are referred to as the second-order displacement fields. 
The imperfection shape is expressed as 

{ u° ) = q°{Ui } + q° q°{Uij} ( i,j = 1, 2, ... ) (2.4.2) 

where q° is the "modal imperfection amplitude" for an imperfection in the shape of {«,}. 

The boundary value problems governing the three families of shape functions are obtained by ex- 
panding the equilibrium equations and boundary condition expressions in terms of the expanded forms 
of the displacements, grouping terms based on their order in the modal amplitudes, and setting to zero 
the imperfections and pressure. The resulting expressions of order 0, 1, and 2 provide the boundary value 
problems governing [uc], { u, } , and {«,,}, respectively. Function { Ul ) satisfies the non-homogeneous 
in-plane boundary conditions, whereas functions {«,} and {«,>} satisfy homogeneous in-plane boundary 
conditions. Each of the three sets of shape functions is discussed in a following section. 
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2.5 Prebuckling Solutions 

, c , A u ™ t 0 1 ° ad " ystem for ^ prebuckling response is prescribed as part of the problem definition 
section 2.3.1). The unit load system consists of some combination of unit global in-plane loads AU 
and Nyct, and unit global length- and width-change measures A it and Av t , consistent with the desired 
form of m-plane loading selected from Tables 1-3. The unit loads are distributed so as to produce uni- 
form normal displacement of both the longitudinal ends and the side edges. It is assumed that the pre- 
buckling response involves only uniform in-plane normal strains in the component plate strips (no plate 
bending). Thus, neglecting rigid-body displacement of the plate strips, the prebuckling solution \u,] 
involves only m-plane displacements: 1 

T 

{ u l) = E U L V L ^3 (2 5 l) 

The nature of the loading and elastic properties of the component plate strips are such that neither shear 
strains not shear stress resultants are induced, so that 


«-{?} 


(2.5.2) 


and 


{AU=[A]{eJ = 


-m 


(2.5.3) 


TTie quantities {ej and {Af t } are uniform over each plate strip, and this unit solution satisfies the equi- 
librium equations (2.1.8). The detailed formulae used for computing the prebuckling solution are given 
in Appendix A. The foimulae are compatible with those used in PASCO [13] for determining the 
prebuckhng response, but additional equations are included for determining the unit global loads N**. 
ana Nyct (required as input to PASCO) corresponding to specified values of A Uc and/or Avz, 

NLPAN, like VIPASA, requires the strains {£<.] but not the displacements {it}. For some doubly 
connected plate-stnp configurations, the method used in PASCO to compute { 14 } results in a violation 
of displacement compatibility between adjoining plate strips. This is because the method only enforces 
compatibility of the longitudinal axial strain e, between the various plate strips. The NLPAN code in- 
cludes an option to enforce displacement compatibility for cases where stiffener flanges are modelled 
as plate strips which are offset from the skin of a panel. In the linear solution, the plate strips are con- 
, ned response. Used in a pre-processor mode, NLPAN provides PASCO input param- 

eters NX(1), NY(1), and FNY(I) [14] which distribute the in-plane loads among the individual plate 
stnps so as to satisfy the displacement compatibility condition. When the prebuckling solution is com- 
puted with the enforcement of displacement compatibility and zero plate bending, there will in general 
be non-zero values for the associated node-line moment resultant contributions Ml, 

2.6 Buckling Eigensolutions 

The buckling equations obtained in the manner described in Section (2.4) are given by 

N v x + N *yi’y + VV W <’»’ = 0 (2.6. In) 

^xy,' x Ny’y + \N Xf V i , xx = 0 

M x x xx + 2jW iy,’x> + M y,-yy + + N y{ W ityy ) = 0 

Where X. is the i* buckling eigenvalue, and functions [A/,] = (X N„ NjY and [MA = [M M M T 
given by ’’ 1,1 J 


(2.6.16) 
(2.6.1c) 
are 


WJ-DOte) 


m =Mk} 
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( 2 . 6 . 2 ) 



where 


{£/} = 




(2.6.3) 


The eigenvalues X, are generally negative, corresponding to a compressive end load on the panel. 

Equation (2.5.1a) differs from the corresponding equation used in the VIPASA analysis 1121 the 
latter equation being 

■*" N x y.,y + = 0 (2.6.4) 

The third term of each of the disputed equations is assumed here to be negligible, based on arguments 
given in Ref. [7], and so equation (2.6.1a) is replaced by 

N *i’* + N *y,'y = 0 ( 2 . 6 . 5 ) 

The generalized side-edge force resultants \f. ) = \j„f y .f u of equation (2.1.10) have contributions 
lA*) \jxufyiMjiU Wk*_\ associated with the buckling eigensoludons, given by 

I y 

w y ‘ 

Hy&Mxyt’x ■*" Myfy ^flyftt’y ) I 

-VM, 

to express the node-line boundary conditions (Section 2.3.3) to complete the 
boundary-value problem specification for the buckling mode shapes. 

"Hie solution of the buckling equations for a linked-plate configuration is performed by the VIPASA 
buckling and vibration analysis code [12]. The r rt eigenfunction has the following assumed form on each 
plate strip (where a phase shift in the x-direction has been applied relative to the conventions of VIPASA 
to provide that w,{0, y) = w,(L, y) = 0 ): 




{«,} = 



£,{y) cos mjixIL 
T| ( <y) sinm,7Cx/L 
<]),<y) sin m.nxJL 



(2.6.7) 


where m, is the integer number of buckling half-waves along the length of the panel for the i* buckling 
eigensolution. When the buckling equations (2.6.1) are expressed in terms of the displacement forms 
given above, the former equations can be reduced to two coupled homogeneous linear second order or- 
dinary differential equations in the functions £,,{y) and T|,(y), and one homogeneous linear fourth order 
ordinary differential equation in the function ^{y). 


For specified values of the longitudinal halfwave number nu, the VIPASA program returns buckling 
eigensolutions in the form of an eigenvalue h and a set of generalized node-line displacement ampli- 
tudes. Using this information along with the general forms of the solutions to the governing ordinary 
differential equations, the functions {^<y)} = [ %<y) n.<y) <t>,<y)] r are obtained in analytic form. Details 
of the procedure are described in Ref. [7] Appendix B. (The following error in [7] is noted. Parameter 
L in equation (BIO) should be given by 
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1 


~k(- x l&+ D * t - D ") 


Furthermore, the term U appearing in the definitions for L and T of equation (BIO) refers to the length 
squared, and not the parameter L defined in the equation above.) For an infinite-length structure sup- 
ported at uniform intervals L along the length, the eigenfunctions satisfy (to the first order) the homo- 
geneous form of the boundary conditions along the boundary node-lines which have been selected from 
Table 2. However, these boundary conditions may not be satisfied for a finite-length structure, in which 
case a correction, described below, is applied to the buckling eigensolution. 


If boundary condition Option 3 for Component 2 in Table 3 (BCVEC(2,IB)=3) is chosen, the cor- 
responding homogeneous boundary condition .is given by V> = 0 and F,, = 0 where Vf and F yi are the 
contributions of the buckling mode to V" and F, of Table 3. The condition F yi — 0 is satisfied automat- 
ically only if the average value F yi is evaluated over an even number of halfwaves (i.e. m, of equation 
(2.6.7) is even). If m, is odd, the condition F yi = 0 may be violated, which indicates that a spurious con- 
tribution to the global load^ component N yC may be present. Thus, a compensating contribution must be 
added. The amplitude of F yi (the first-order contribution to the load F y appearing in Table 2) is deter- 
mined from the mean value of the stress resultant N yi in the panel skin along one of the two boundary 
node-lines: 


^ = »A.U 


where 


N y _\ m =-j- rVi 

y ‘ "i L J 0 y ' »i 


dx 


( 2 . 6 . 8 ) 


(2.6.9) 


and rti designates a boundary node-line location. The global load contributions Nya is then defined as 


NyG =N y I 

y^i it 


( 2 . 6 . 10 ) 


The associated load component in the x-normal direction is denoted N Mi , and this component is zero by 
virtue of the harmonic form of N„ (/V* = 0 at x = 0, L): 




( 2 . 6 . 11 ) 


= 0 


A modified buckling eigenfunction { u , } is defined, having the form 

{«,) = {«,} +^,{«b) (2.6.12) 

where the shape function {ue} has the same general characteristics as the unit linear unbuckled solution 
{«£.} (see Section 2.5), except that the former corresponds to a unit load system 

N yG = NyG B , N xG = N xGb = 0 (2.6.13) 

Modified global load contributions, N^, and Mc„ are defined, corresponding to the modified 
eigenfunction {«,). In view of equations (2.6.11) and (2.6.13), 

N xG= 0 (2.6.14) 

Parameter B , is chosen so that the following homogeneous boundary condition is satisfied: 
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(2.6.15) 


"yorNya+B^ 

= 0 

with the result is that 

B i = - N yc/Nyc B (2.6.16) 

For modes which have non-zero values B„ the eigenvalue A, associated with {«,} will be slightly different 
from the original eigenvalue A,. NLP AN computes and prints the values A,. 

If the value of BCVEC(2,IB) is 3 for one boundary and 2 for the other (see Table 3), then N& may 
be non-zero on one boundary and zero on the other. This discrepancy is the result of the nature of 
VIP ASA buckling mode shapes. In NLPAN, the modified eigenfunction { u, } is used only if 
BCVEC(2,IB)=3 for I B= 1 and IB=2 (both y-normal boundaries). 

2.7 Second-Order Displacement Fields 

The equations governing the second-order displacement fields, obtained in the manner described in 


Section (2.4), are given by 


+ ^xy t j’y + ~2 ■*" = 0 

(2.7.1a) 

Njcy-'X 4" Ny^yy + ij'XX + = 0 

(2.7.1ft) 

^XtfXX + 2Afx y..,xy + Myjyy + X^IJ^XX + ^y^ijiyy) 

+ ~2~ t + ^xy^yy) ^ (fixyWjt x + ^y^j f y)'y 

(2.7.1c) 


+ (N x yVfi x 4" ^xy^i’y) 'x 4~ (fix y^i'X ^ — 0 

where (V y ) = N yiJ and { Af t> } = M }iJ are given by 

i N ij) = U] (Cy) , [My } = [D] { k ij] (2.7.2) 

where 

ij'x 4* w i<x w j*x) f w ij<xx ^ 

v ijty -5{Ui,yUj,y + Wi, y Wj,y) >■ , (lCy}=— ^ij'yy J (2.7.3) 

U ij,y + Vy.j. + ,5{\Vl, x Wj,y + yVj,yWj,^) J ^ 7.W J 

The in-plane load parameter A has been replaced in the above equations by a fixed reference value A*, 
consistent with the desire to obtain a load-independent set of shape functions [lUj). In Section 4, it is 
argued that the load-dependent terms can be eliminated from the above equations, so in the implemen- 
tation of NLPAN, A* is set to zero. (The term A which naturally arises in deriving equation 
(2.7.1a) has been omined, consistent with the deletion of the related term from equation (2.6.1a).) 

Equations (2.7.1) involve known eigenfunctions {«,) and {«,}, and unknown function [u,,). Let m, 
and m, denote the number of buckling halfwaves in the ^-direction for eigenfunctions {«,•} and (uj, re- 
spectively. Separation of variables can be used to convert the trio of partial differential equations (2.7.1) 
into two trios of ordinary differential equations in the variable y, by assuming the following functional 
form for {%} [15]: 
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(2.7.4) 



where 


{ m i + rrij a = 1 

m, - rrij a = 2 


(2.7.5) 


The transverse displacements v v and w v do not satisfy the same boundary conditions as v, and w, at 
x = 0,L. This issue is discussed in detail in Section 4. 


Using the assumed forms for {«,,} and {«,>} along with trigonometric identities, the partial differ- 
entia] equations (2.7.1) are converted into two uncoupled trios of non-homogeneous ordinary differential 
equations, one trio governing {£,«} = . the second trio governing {^,} = ft* q*, fcj . 

Both equation trios have the following general form, where the subscripts a and ij are omitted: 


+ C& + C 3 r|' = F(y) 

(2.7.6a) 

Dg + D 2 T]" + D 3 T\ = G(y) 

(2.7.6b) 

£,<{>"" + £^<t>" + £^<|> = //(y) 

(2.7.6c) 


where the subscripted C’s, D's and Fs are constant coefficients, and the iunctions F(y), G(y), and H(y) 
have a quadratic dependence on the components of functions {^(y)} and (^(y)} and their derivatives. 
For general configurations, all three equations are coupled through the boundary conditions along 
node-lines where non-coplanar plate strips join. The expressions for the coefficients and nonhomogene- 
ous terms in equations (2.7.6) are presented in Section 3.5 of Ref. [7]. 

The homogeneous form of the boundary conditions discussed in Section 2.3 are applied to the 
second-order displacements and the associated generalized node-line force resultants. These boundary 
conditions complete the specification of the boundary-value problem governing the functions {^a,Xy)}. 
Detailed boundary condition expressions are presented in Section 3.5 of Ref. [7], Because of the com- 
plexity of the expressions found in IX^y) G aj (y) //*,<y)] of equations (2.7.6), a finite difference analysis 
is employed to obtain values of the functions {^{y)} at a finite number of evenly spaced points across 
the y-domain of each plate strip. The finite-difference solution procedure is presented in detail in Ap- 
pendix C of Ref. [7]. 

The functions { 14 ,} satisfy the homogeneous form of the conditions along the boundary node-lines 
which have been selected from Table 3. At the ends of the structure (x = 0, L), {«,,} satisfies the homo- 
geneous form of the displacement-control (CONTRL=’D' in Table 1) boundary conditions. If load 
control (CONTRL='L' in Table 1) is used, modified second-order displacement fields {«,,} are computed 
which satisfy homogeneous force boundary conditions at the longitudinal ends. The modified function 
has the form 


{“<>} = («,)} +Aj<{u a } ~ {«/.)) (2.7.7) 

where {14} = v L 0 ] is the unit prebuckling solution discussed earlier, { 4 ,} = [0 v A 0 ] T is a unit 
unbuckled solution of the same nature as { 14 } except that it corresponds to a unit global load NyCA where 
the corresponding longitudinal strain e*,, is held to zero throughout the structure, and A„ is a constant 
coefficient to be determined. The source and significance of the modifying terms are discussed below. 

The second-order displacement fields give rise to associated second-order contributions to the global 
loads Njg and N^, according to the relationships 
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(2.7.8) 



(2.7.9) 


where the quantity inside the parentheses in equations (2.7.9) is the value of N# in the skin at one of the 
boundary node-lines. For the case of bi-axial loading (Load Case A, Table 2), the homogeneous form 
of the boundary condition which is satisfied is that of Option 2 or 3 for Component 2 in Table 3, so that 

(2.7.10) 

The modified field {«,>} must contribute nothing to either N* or A/*, so the following conditions are 
imposed: 


=0 

A^ = 0 (2.7.11) 

With substitutions, the above two equations take the form 

N xGij + A ij( N xG A ~ NxG l ) = 0 

Aij(NyG A - NyG,) - 0 (2.7. 1 2) 

The value of N& is selected such that 

N yG A = NyG L (2.7.13) 

and A,j is then given by 

A u = N xgJ(H x g l -N xGa ) (2.7.14) 

For load-control cases in which the desired node-line boundary condition given by Option 1 for 
Component 2 in Table 3, the condition V* = 0 along the boundary node-lines is automatically satisfied, 
but it is still necessary to use the modified function {jjy} satisfying 

N *u = 0 (2.7.15) 

This is accomplished by using the expressions for {«,,} and Ay given in equations (2.7.7) and (2 7 14) 
respectively, but setting to zero the quantities {u x }, A^, and N^. 

2.8 Stationary Total Potential Energy Formulation 

2.8.1 Expansion of the strains, curvatures, and stress resultants. The form of the displacements 
presented in equation (2.4.1) is restated with the modified functions {«,} and {«,,} of equations (2.6.12) 
and (2.7.7), respectively, replacing the functions {«,} and («,>): 

{«} = Mu L ] + qdii) + MjiHij) ij= 1, 2, 3,... (2.8.1) 

where 


(«,} = K) +s t [u B } 

(«//) = (««>} +Ajj{{u A } - {u L }) 
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( 2 . 8 . 2 ) 



The generalized in-plane load parameter X controls either end displacement or end load, consistent with 
the selection of parameter CONTRL from Table 1. Simplifications to equations (2.8.2) which arise in 
some situations are summarized in Table 4. 


The mid-surface strains and curvatures of equation (2.1.3) are expanded by using the expanded form 
of the displacements given in equations (2.8. 1-2), giving rise to the following forms: 

{e m } = X(e L } + (<7, - <7f){e,} + (q^ - q°qf)[i t j) 

+ (< 7 .< 7 / 7 * - (£,>*} + (qfl/lkQi ~ 

tO =(?,- <7 D { *,} + (q - q?q °) ( *<>} (2.8.4) 


where 


{£.} = (£,) +5,(6^} 

{£«>} = {£«>)+ A><{eyt) -{££.}) 


(2.8.5) 


The terms {e t }, {e,}, {£,,}, {k,}, and {k,;} appearing in equations (2.8.3)-(2.8.5) have been presented in 
terms of displacement functions in equations (2.5.2), (2.6.3), and (2.7.3). The remaining terms are de- 
fined by 

r 0 1 

le B }= < v Bfy J- {£,!}=< v*,, V 

^ ^ J (2 8 6) 

f v bx v jbx + w i>x™jhx 1 f ■5(v i j, x v kJ , x +w i j, x w kl , x ) "I 

— U i<y U jfry W i’y W jk’y S’ {£(/*/) = % ■5( u ij<y U kl'y + W ij’y W kl’y) S 

V W nx W jk'y + w i>y w jbx J ■^( w ij’x W kl'y + w ij'y W kJ'x) J 


Stress resultants { A^} and [M] are linearly related to the mid-surface mechanical strains and cur- 
vatures, respectively, through equations (2.1.11), and thus have the same expanded forms: 

{IV} = X{N l } + {q i - q°)[N t ) + (qq - q°q°){N tJ ) 

+ (qaflk - q°qjq 0 k){^ijk) + - qUjqtqt)^^) u ‘ 

\M) = (q, - q°) {A/,} + (q,qj - qfqfr[Mg) (2.8.8) 


2.8.2 Edge-load contributions to the total potential energy. Global load and displacement quantities 
are used to express the contributions of the in-plane loads to the total potential energy. Consider the 
in-plane displacements, which have the expanded form 




Ui+BiUg 
Vi + BiVg 



Uij + A i} { 0 - U[) 
v y + A tji v A - vj 


} 


(2.8.9) 


At the longitudinal ends of the structure, «,y is identically zero. It is also assumed that non-zero values 
of the axial buckling displacements w, = w, + B,Ue at the ends of the structure are due to rotation of the 
ends about the generalized neutral axis without any associated change in length. The effective change 
in panel length, Aw, is then 


Au = u 


x = L 


- u 


x = 0 


= (k- qfljA L] )&u L 


( 2 . 8 . 10 ) 


The mean change in panel width, Av, is given by 
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( 2 . 8 . 11 ) 



- XAv l + Av,- + B; Avg) + <?tf y [Av,y + A ij {^v A - Avj] 

where over-bars denote the mean quantities taken over the length of the panel, and V' and V* are the 
displacements of the two boundary node-lines in the global y -direction. 


The global in-plane load components N& and N# are defined by 



( 2 . 8 . 12 ) 


(2.8.13) 


where «i is a boundary node-line, and N, is measured in the panel skin (if present) at the boundary 
node-line. The expanded form for {N} of equation (2.8.7) is inserted into the above two equations to 
obtain the following expanded form for the global in-plane loads: 




+ (?, +q°) 


{ 


0 








} 

} 


+ 


f OCq, 4 ) 
1 <Xq?) 


} 


(2.8.14) 


The higher order contributions inside the last set of braces are neglected. 


For the case where load control is used for the in-plane loads (CONTRL='L' in Table 1) the con- 
tributions of the in-plane loads to the total potential energy (see equation (2.2.3)) are now expressed in 
terms of the global parameters as 

a a 

-BN^Au - Nyc Av dy (2.8.15) 

J 0 


Simplifications apply to the expanded forms of the terms appearing in the above expression because load 
control is in effect. Load component N& becomes 

A 

n xC = ^xG l (2.8.16) 

so that 


-BN xG Au = (-X 2 + qfljkAipEt (2.8.17) 

where 

Et = BN xGL Au L 

For Option 2 of Table 1 (CONTRL=’L') the load N# is given by 


(2.8.18) 



- ji WyG Av dy = -X£>[ - q,X(Dt + Bp L B ) - q,q^ + - £>[)] 


( 2 . 8 . 20 ) 


where 

Dt = LNyc L Av t 
D b =L NyQ L Avg 

Da = L NyQ L Av^ 
= L Nyc L Av, 
D^LN^j 

For Option 1 of Table 
contributes nothing. 


( 2 . 8 . 21 ) 


(CONTRL=’L') t Av is zero, so that the second term of expression (2.8.15) 


2.5. J Expansion of the total potential energy. The total potential energy expression, equation (2.2.3), 
is rewritten here, reflecting the use of the global load and displacement parameters to represent the 
contributions due to applied edge loads: 


p 

* = X(7 y({N}V} + {A/} T {K c }>i4- 

p=\ Ka J a / p 

-P 

J o 


— B N^q A u — NyQ Av dy 


( 2 . 8 . 22 ) 


The pressure load Qp acting on plate strip p is proportional to a global pressure load parameter p, such 


Qp~$Qp (2.8.23) 

A 

where Q, is a specified unit value for The total potential energy is now evaluated in terms of the 
expanded forms of the displacements (equation (2.8.1) ), the mid-surface mechanical strains and curva- 
tures (equations (2. 8.3,4) ) and the stress resultants (equations (2. 8. 7, 8) ) and the substitutions of 
equations (2.8.17) and (2.8.20) are made. The following equation is obtained: 

n = Constant 

+ { - PC? - qp - qjqpjt - q?qtqfcj u ) 

+ *T ti- I&v - 

"*■ + Cij — q t C ijk ) 

+ + Cl jk + -i- C %) } + 0(< ?1 5 ) 


where the "Constant" terms are those which do not depend on the modal amplitudes. The coefficients 
with over-bars are given by 


£? = (Ct + B,ck) - {D[ + Bpjf) 

+ B,C f + BjCf + Bpfi 

% = ICti * - Ct)l - 1 Dfj + A.,{D l a - (Dt + £?))] 

Cij = Cij + - C k ) + B k C t j + AijB^cfl - C B ) 

cjjk ~ Cp + B[Cf jk 

c“=c“j-A u dt- a A + W -icji+cb 
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The unbarred coefficients are "primitive" coefficients. Coefficient is defined in equation (2.8.18), 
til, Db, Da, Df, and Df, are defined in equations (2.8.21), the pressure-related coefficients are defined by 



and the remaining primitive coefficients are defined as illustrated here by two examples: 


P 

C* = S(J + W k ) T [<ij\)dA 

p=l' A 

The primitive coefficient Cj which arises in determining equation (2.8.24) is omitted because it is 
identically zero. Contributions to n which are of order five or greater in the modal amplitudes are 
henceforth neglected. 

The displacement shape functions and their derivatives are known analytically in the x-direction, 
and are known at discreet points across the y-domain of each plate strip. Thus, in evaluating the various 
integrals which define the primitive coefficients, analytical integration is performed in the x-direction, 
and numerical integration (using Simpson’s rule) is performed in the y-direction. Numerical integration 
is performed using quantities evaluated at the discrete points in the transverse plane used for computing 
the finite-difference solution for the second-order displacement field. 

A stationary total potential energy condition is imposed in order to obtain expressions governing 
equilibrium: 

5 * = 8 *'(l?“) (2.8.28) 

= 0 

Each of the "virtual modal amplitudes," Sq, (i = 1, 2, ...), appearing in equation (2.8.28), is both arbitrary 
and independent, and thus, in order to satisfy the equation, each expression dn/dqt must independently 
be zero. With the selection of a finite basis of M buckling mode shapes for use in the analysis, the 
following set of simultaneous algebraic equations is obtained: 


). 


(2.8.27) 


j*- = (*# + PC? - q° C\ - qjq° k cf - q'qWC? 1 ) 

+ qfljfy + pC? + C, y - <7*C* - qUiC“) 

+ q J q k (\C^ + C ljk -q?c! jk ) 

+ + c lJkI ) 

= 0 0 = 1,2 M) 

where the modified coefficients which appear are given by 


(2.8.29) 
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C{ k 

c{ u 

c*= 

z 

A 

0/A : 


= c 


= c; 


■ju 

2C k 

2C% 

■■ICt + cik 


c!j k =c! Jk +2^ jik 

Ci jkl 

C ijkl 


= 2(C£* + Cf w ) 

= Ctjk + 2Cy i * + Cju + C,* 7 + ^ 


C^ = 2Cj 

cf=cf 

3 = 3 

C[ jk = c[ ]k + 2C^ 

cF=-cp 

cf = -2Cf 


(2.8.30) 


Equations (2 1.8.29) are nonlinear in the modal amplitudes q, (j = 1,2, ... , A/) and linear in the two 
load parameters A, and p. The algebraic equations have been derived for general values of the modal 

imperfection amplitudes q? (i = 1,2 Af), and thus the equations can be used to explore the structural 

response in the presence of a variety of imperfection shapes and amplitudes. 


2.9 Solution of the Nonlinear Algebraic Equations 

2.9.1 Normalization of the variable parameters. Before applying solutions strategies to the set of 
nonlinear algebraic equations (2.8.29), the variable parameters are normalized so that they take on values 
of order of magnitude unity. The parameters to be normalized include the set of modal amplitudes { q) 
and the in-plane and pressure load parameters, X and p, respectively. Each of these is discussed below! 

The buckling modes {u,} are normalized before their subsequent use in other parts of the NLP AN 
anaiysK, such that the maximum displacement for each mode anywhere in the structure is equal to a 
specified reference thickness /u/. Thus, the modal amplitudes appearing in equations (2 8 29) are con- 
sidered to be normalized. 


load parameter X is normalized by its critical value associated with primary buckling, 
Ki. Therefore, the normalized in-plane load parameter X is given by 


X = X/Xi 


(2.9.1) 


The pressure load parameter p is normalized by a reference value P„ 7 so that the normalized pressure 
load parameter p is given by 


P = fMW 


(2.9.2) 


The reference value is selected to be the value of p which produces a deflection response, in a single- 
mode NLPAN analysis, having a maximum amplitude in the structure equal to the reference thickness 
Thus ’ in a sin & le mode analysis of the perfect structure using mode number r, if p is set to p r>/ and 
A, is set to zero, then the modal response is \q r \ = 1. Applying equations (2.8.29), it can be determined 


fW~ 



( C rr + C rrr + C rrrr ) 


(2.9.3) 


where the coefficients shown are those defined in equations (2.8.30) with the hats omitted, and where 
cretticient C% which would generally appear has been assumed to be negligible compared to coefficient 

c ; • . number r is selected to correspond to the mode which, among those in use, has the co- 

efficient Cf of largest magnitude. 


Equations (2.8.29) can now be expressed in terms of the normalized parametere X and p The re- 
sulting equations are identical in form to equations (2.8.29) except that X and p replace X and p, and the 
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constant coefficients have been transformed in the obvious manner to account for the definitions of A 
and p. 


2.9.2 Specification of shape imperfections. After the normalization of the variable parameters is 
completed, the general form of the nonlinear algebraic equations is unchanged from the original form 
given in equations (2.8.29). To simplify notation, the notation of equations (2.8.29) is used in this sec- 
tion, where it is understood that the parameters and coefficients which appear are the normalized ones. 
Once the values for the modal imperfection amplitudes, [q°], are specified, equations (2.8.29) can be 
written as 

(Cf + xcf- + pc?) + qff*. + AC? + pc?) 

+ <iM C ijk + ^iik> + 9fh^iC tjkl + AC?*,) = 0 (i = 1 , 2, ... , M) (2 ' 9 ' 4) 


where 


c? = - - - qf tf cf - tffcVcf 

c°jk = c ijk - q° c\jk 


(2.9.5) 


and where the hat symbol 0), appearing over the coefficients in equation (2.8.29), has been dropped here. 

Given the shape-imperfection field (r) the effective modal imperfection amplitudes can be deter- 
mined by making use of the orthogonality, condition for the buckling mode shapes. With the i* buckling 
mode shape denoted by {u,} = []w, v, w,J , this orthogonality condition is given by (see Appendix B 
equation (BIO)) 


<{«;}, £({«,})>= | 


i=j 

i*j 


(2.9.6) 


where the inner product on the left-hand side of the above equation is given by 


<{“>}. m «;})>= 



"I" kVj{N x \V^, zx 4" NyW^yy 




(2.9.7) 


where p is the plate strip index number, P is the total number of plate strips, A is the planform area of 
a plate strip, and and are the in-plane stress resultants corresponding to the unit pre-buckling 
solution. Express {u°} as a series in the buckling mode shapes, 

M 

1“°) (2-9.8) 

1=1 


and evaluate the following inner product with consideration of the orthogonality condition: 

M 

< { u°) , L( { u, :} )> = { Uj] , L( { «,))> 

j= i 

= q°bi (i = l,2 M) 

where no summation over i is intended. Hence, the modal imperfection amplitudes are given by 


(2.9.9) 


q° = ±<{u 0 m[ Ui ))> 0 = 1,2 M) 


(2.9.10) 


The computation of the coefficient b, is performed in the NLPAN program for all modes. The evaluation 
of the inner product of equation (2.9.10) is more challenging, since the shape-imperfection field is not 
generally known as a continuous function, but is more likely known as a set of values at discreet points, 
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or is known along a set of discreet lines on the surface of the structure. No attempt has been made in 
NLP AN to automate the computation of qt based on a measured imperfection field. 

2.9.3 Asynchronous control of in-plane and pressure loads. For the general case of combined in- 
plane and pressure loading, the nonlinear algebraic equations governing equilibrium (equations (2.9.4)) 
contain the two load parameters X and (3. A single load parameter. A, is used when applying the non- 
linear solution strategies discussed in the next section. This section describes a procedure for relating X 
and P to A in a way which permits the modelling of asynchronous application of the two types of 
loading. 

A series of load ranges is specified in terms of target values for X and (3: (0, 0), (Xi, PO, (Xz, P 2 ), 
... . Over the load range, X and p vary linearly with A as A varies from 0 to 1 : 



0< A< 1 


(2.9.11) 


Equation (2.9.1 1) is used to express X and P in equation (2.9.4) in terms of A, and terms are regrouped 
to obtain equations with the following form which governing equilibrium over the k* load range: 

(Q + a ch + qfCij + AC,”) + q flk {C ljk + ACj*) 

~ (2 9 12) 

+ 4 'W( C W + = 0 (i = 1 , 2, ... , M) 

Details of this procedure are given in Appendix D Section D.l, where it is noted that there are some 
obvious differences in notation. Equations (D3), (D5), and (D6) of Appendix D correspond to equations 
(2.9.4), (2.9.11), and (2.9.12), respectively. 

2.9.4 Nonlinear solution strategies. The solution of the system of nonlinear algebraic equations 
governing equilibrium (equations (2.9.4) ) can in many cases be performed using Newton-Raphson it- 
eration. In some important cases, however, limit-point behavior or solution branching is encountered, 
and Newton-Raphson iteration breaks down in the vicinity of associated critical-stability points. Conse- 
quently, advanced solutions strategies for the analysis of nonlinear equilibrium behavior have been im- 
plemented in NLPAN. The implementation of the strategies is the subject of Appendix D, which was 
extracted from [16]. The two basic strategies incorporated are the arc-length control method popularized 
by Riks [17] and the equivalence transformation method developed by Thurston [18]. An integrated 
procedure is obtained which allows equilibrium solution paths to be followed past limit points, and 
through regions of complex solution branching involving modal interaction. 

The implementation of the solution strategies has been found to be fairly robust, with a few quali- 
fications mentioned here. The solution strategies require the specification of several parameter values 
which, for example, control the solution step size, and determine cutoff values used for determining the 
proximity and classification of critical stability points. On occasion, the solution strategies are "fooled" 
into making an improper diagnosis, and the solution stepping will proceed along an inappropriate sol- 
ution path. An adjustment of the input parameter values will generally fix this problem. This type of 
difficulty arises when very small modal imperfection amplitudes are used, or when there are very small 
effective imperfection values arising due to numerical error. The use of significant modal imperfection 
amplitudes (q° > .01) tends to improve the reliability of the solution procedures. 

It has also been found on occasion that a solution path will lead into a web of unstable equilibrium 
paths from which the solution strategies are unable to guide the analysis. This behavior occurs in the 
highly nonlinear regime (such as deep postbuckling), and is believed to be associated with the use of 
an insufficient set of buckling modes as shape functions. If the tangent stiffness matrix becomes singular 
for the second eigenvalue (indicating that this troublesome behavior is being encountered), the NLPAN 
analysis is terminated, and a warning message is issued. 


20 



3. EXTENSIONS TO THE METHOD OF ANALYSIS 


3.1 Dynamic Snap Analysis 

The nonlinear analysis of NLPAN is designed to predict the equilibrium response of a structure 
which is subjected to a varying load level, and thus there is an assumption made that the time rate of 
loading is sufficiently small so that inertial effects are negligible. The methods of Riks [17] and Thurston 
[18] have been implemented for following the equilibrium solution paths past limit points and bifurcation 
points, and the ability to track unstable equilibrium paths has been achieved (see Appendix D). While 
a physical structure will not follow an unstable equilibrium path, following such a path with the analysis 
will sometimes reveal a second path or path segment exhibiting stable equilibrium. In this situation, the 
physical structure is assumed to exhibit a dynamic snap from the original stable path to the second stable 
path, departing the original path at the point where the equilibrium becomes unstable. 

In some test cases modelled using NLPAN, the analysis follows an unstable equilibrium path to 
secondary limit points or bifurcation points where there is no transition back to solutions of stable 
equilibrium. This situation is detected when the second eigenvalue of the tangent stiffness matrix goes 
to zero. The process of trying to navigate through such a complex regime of solution paths (in search 
of the desired path of stable equilibrium) is difficult to automate, and is also of questionable wisdom 
since the solution paths represent structural behavior which will never be encountered. The snap phe- 
nomenon is a dynamic event, and by analyzing it as such, the difficulties encountered in following 
equilibrium paths are avoided. This reasoning provided the motivation for incorporating a dynamic 
analysis capability along with the static analysis capabilities already in place. 

In the time since this undertaking was proposed and initiated, a change in thinking has taken place 
in the mind of this investigator. The primary mode of response for which the dynamic solution capability 
was desired is that of the change in waveform which occurs in structures dominated by local- buckling 
behavior (as opposed to global column- or wide-column-type buckling). A study of the literature has 
shown that the accurate representation of secondary instability behavior generally requires the use of a 
relatively large number of appropriately selected mode shapes as global shape functions (see, for ex- 
ample, [19]). While NLPAN employs rigorous criteria for determining critical stability points, it in- 
corporates no strategy for selecting the crucial modes required for accuracy, and NLPAN is also limited 
to a relatively small number of included modes. For these reasons, a dynamic analysis capability has 
been judged to have relatively little practical value in the NLPAN code. The dynamic analysis capability 
detailed in this section has been coded in NLPAN to the extent of computing the generalized mass and 
damping coefficients, but the full solution procedure has not been implemented. 

The theory for the dynamic snap analysis is presented here. The goal of the dynamic analysis is 
to locate the new equilibrium position sought by the structure during the snap, so the accuracy of the 
dynamic analysis is not considered to be of extreme importance, and some simplifying assumptions are 
thus made. The dynamic analysis is perfoimed at a fixed value of the generalized load parameter A 
(defined in Section 2.9.3), namely the value at which the equilibrium path becomes unstable. Motion- 
damping forces are incorporated. In the following sections, the equations of motion are first developed, 
and then a solution procedure is described. 

3.1.1 Equations of motion. Hamilton's principle (as it is presented in [20]) is applied to obtain the 
equations of motion. Let T be the kinetic energy of the system, and let SW be the external virtual work 
done on the system, which is broken down as follows: 

8W = -5ti + 8W„ c (311) 

where 5 tc is the first variation of the total potential energy for the conservative elastic system and the 
specified loads (see equation (2.8.28)), and 6VV* is the external virtual work of non-conservative forces, 
which in this situation are the motion-damping forces. Hamilton’s principle can then be stated as 
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c ' 2 

J (dT+8W)dt = 0 subject to (8m) =0 at t= ^ 


where /, and h are two reference times, t,<t 2 . Using equations (3.1.1) and (3.1.2), Hamilton's principle 
is re-expressed in its negative form as 


J (5tc - - 87) dr = 0 subject to (6 m) =0 at / = /,, ^ 


(3.1.3) 


The kinetic energy is taken to be that associated with motion of the mid-surfaces of the plate strips 
in the transverse (y- and z-) directions. This is written as 


■Mi 


m(v 2 + vv 2 ) dA 


(3.1.4) 


where m is the mass per unit area of a plate strip, p is the index number of a plate strip, P is the total 
number of plate strips, and / denotes the derivative of function /with respect to time. Using integration 
by parts and imposing the constraints {6u} = 0 at t = f, and t = h , the following expression is obtained: 


J (57) = - ^ ^ J m(v8v + w5w) dA^ dt 


(3.1.5) 


In evaluating the right-hand side of the above equation, only the first-order contributions to the dis- 
placements v and w (see equation (2.4.1)) are taken into account, to arrive at 


% 

✓ 

Qj = £ ^|m(v,v, + w ’-Wj) dA^ 


(3.1.6) 


(3.1.7) 


A non-conservative motion damping force is assumed to act normal to a plate surface in the di- 
rection opposite its normal velocity. The force per unit area is assumed to equal the normal velocity, 
w, times a damping coefficient, p, which has the units of force per unit area per unit velocity. The 
non-conservative virtual work is thus given by 


-la 


]\w8wdA 


(3.1.8) 


In evaluating the right-hand side of the above equation, only the first-order contributions to the dis- 
placement w are taken into account, to arrive at 


M M 


(3.1.9) 


<=i j = l 


where 




(3.1.10) 


22 



With the generalized load A held constant at the critical value A* (the value at which equilibrium 
becomes unstable), bn can be written as 

Stc = + qf tj + qfl k C ijk + (3.1.11) 

where summation over repeated indices is implied, and the coefficients appearing in equation (3.1.11) 
are related to the coefficients defined in equations (2.9.12) by 

C, = (A*cf + C°) C, Jk = (A* c[ !k + C° k ) 

C s =<A*C$ + (~) c iju = (A* cjju + C?j U ) 0U2> 

Equations (3.1.6), (3.1.9), and (3.1.11) are now substituted into equation (3.1.3): 

J & 7 ,<C/ + qfij + qfltCijk. + <lflkRfiju + qfij + qfij) dt = 0 , subject to { 6w } = 0 at t = f, . ^ (3.1.13) 

h 

The time interval [ft, fc] of equation (3.1.13) is arbitrary, and each of the virtual modal amplitudes 
oqi (i =1,2, ...) is both arbitrary and independent, so that the equations governing motion at a fixed level 
of the generalized load can be reduced to 


qfy + qf^ + qfj + qfl k C ijk + qfl^f ijkJ + C, = 0 

0 =1, 2, ... , M) < 3 - U4 > 

Initial values of the generalized coordinates and their time derivatives, [q),~ 0 and {^}, =0t respectively, 
must accompany the equations of motion. One approach is to set the former equal to the critical stability 
solution {< 7 *} and set the later to 

{<7}/=o = e(<t> 1 } (3.1.15) 

where is the eigenvector corresponding to the zero eigenvalue of the tangent stiffness matrix, and 
e is a small non-zero number. 


3.1.2 Solution of the equations of motion. The Newmark direct integration procedure, as presented 
in Ref. [21], is used to obtain a series of discreet solutions to the initial value problem stated above. 
A uniform time increment. At, is used. Index n denotes the current time step number, at which the sol- 
ution is known, and n + 1 is the index number of the following time step, at which a solution is sought. 
The Newmark procedure uses two somewhat arbitrary parameters, a and p, that determine the exact form 
of the time-derivative approximations used. The following parameter values are used here: a = 1/2, and 
P = 1/4. This leads to the following time-derivative approximations: 



(3.1.16) 

«r'=<i?+Attf+f (£+£+') 

(3.1.17) 

Equation (3.1.17) is solved for 


"f* + 1 ^ n + 1 n ■ n ..« 

<7; = a oQt -a Q q t -a x q { - q t 

(3.1.18) 

where 


4 4 

a ° = At 2 * 1= A 7 

(3.1.19) 

By expressing equation (3.1.14) at time step number n+\ and eliminating [q" 
equations (3.1.16) and (3.1.18), the following equations are obtained: 

‘} and + using 
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(3.1.20) 


(a 0 C tJ + a&Cij + C^qJ + 1 + * l q n k + 1 

0 =1,2, ... ,M) 


. r* n + 1 n + 1 n + l 

+ <7* <7/ 


= Ei 


where 


£» = - Q + Q ,fq 0 qj + a x q* + qj) + C^aJ^q- + qj) 


(3.1.21) 


and 


*h = ^~ (3.1.22) 

Equation (3.1.20) serves as the basis for determining the solution {< 7 " + l ) in terms of known quantities. 
Once {< 7 "*'} is known, the time-derivatives {< 7 * +1 } and can be determined by applying equations 

(3.1.18) and (3.1.16) in sequence. 

The solution {< 7 "*’} can be determined using Newton-Raphson iteration. Let {q r ) be the r* estimate 
for {< 7 ** '} in the iterative solution procedure. The residual error vector and the tangent stiffness matrix 
for the r* solution, denoted {/?”}, and [AP], respectively, are given by 

Ri — (aoQj + QcfoCij + Cy)qj + C yktfj q^ + CijkfijQkQt ~ £,■ (3. 1.23) 


and 



{/} 


a <fcij + ■*" Qj + qkiPijk ■*" Cjy) ^ <7**7/ (f' ijkl ^ ^ikjl + ('i kip 


(3.1.24) 


The improved solution, is then given by 

W r+1 } = W) + [&q) (3.1.25) 

where the correction, {A*?}, is given by 

{Atf}=-[r:fV} (3.1.26) 

Iteration is continued until the correction { A^} becomes negligible compared to the solution {^ +1 }- 

The damping coefficient pi and the time increment At should be selected based on characteristic 
aspects of the dynamic response. Criteria for selecting these two parameters are not considered here. 


32 Thermal Loading 

The capability of modelling thermal loading has been incorporated in NLPAN for the case where 
each plate strip is subjected to a uniform temperature. ("Temperature" as used here refers to a change 
in temperature relative to a reference value). Each plate strip may have its own temperature value, but 
the temperatures in all plate strips are proportional to a single thermal-load parameter. Thermal loading 
may be used in conjunction with any option for control of the generalized in-plane loading (see Table 
1). Thermal loading is controlled independently of the generalized in-plane loading. 

This section is devoted to describing an additional contribution to the expanded form of the dis- 
placements (equation (2.4.1) ) necessary for simulating thermal loading. The incorporation of the ther- 
mal load contribution into the final problem formulation is described in Section 3.5. 
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A unit temperature T is specified for each plate strip, and this value may be different from plate 
strip to plate strip. The thermal load parameter y is used in the nonlinear analysis, and the temperature 
change experienced by a given plate strip is then T=yT. 

When non-zero thermal loading is present, the expanded form of the displacements on each plate 
strip has an additional contribution y{u r } as shown here: 

{«) = M^l) + Y(«r} + q ,{«,} + qfl Jj{Uij) ( 3 . 2 . 1 ) 

where {«,-} = [ur v,-0] r is a displacement solution corresponding to the unit temperature system. Sol- 
ution [u T ) is a linear, unbuckled solution which satisfies the homogeneous in-plane boundary conditions. 
As with the unit solution {wt} described in Section 2.5, the bending in the component plate strips is 
constrained to be zero for the solution {« r }. 

The expanded form of the mid-surface strains has a new contribution y{e r } shown in the following 
equation: 

{e} =X{e L } +y{£t-} + (<?, -<??){£,} + ... ( 3 . 2 . 2 ) 

where 


{ £ t) - 



(3.2.3) 


Equation (3.2.2) gives the mid-surface strains relative to the imperfect reference configuration. For a 
given plate strip, the in-plane mechanical strains associated with the unit response, {£*}, are given by 
the difference between the actual strains and the thermal strains for an unconstrained plate strip: 


ter) — ter) - (3.2.4) 

where {a} = [ 0 * 0 , 0 ] are the coefficients of thermal expansion for the plate strip. Therefore mid- 
surface mechanical strains {£"}, which compensate for both imperfections and thermal loading, contain 
the contribution y{£?): 


(£ m ) -\{Ej +Y{£t} +(<?, -<??){£,■} + ... (3.2.5) 

The in-plane stress resultants reflect thermal loading by the presence of a new contribution, included in 
the following equation: 


(A) = X[N l ) + y [N t ] + % - + ... 


where 


(3.2.6) 


{A/ r } _ [ 4 ] {e™} (3.2.7) 

and [/l] is defined in equation (2.1.12). The formulae for determining the complete unit thermal sol- 
ution are developed in detail in Appendix A, Section A.2. 


3J Rotationally Elastic End Support 


Rotationally elastic support of the longitudinal ends is simulated using linear rotational springs lo- 
cated at discrete points at the ends of the panel. Spring locations are specified in terms of points on the 
cross- section, denoted as yt, where k is the index number of the spring, and y is used here as a gener- 
alized transverse in-plane coordinate (ie. y signifies both a particular plate strip and a point along its local 
y-axis). Springs of specified strength are then simulated at these points on the cross-section of both 
longitudinal ends. A spring with rotational stiffness K m resists duldy (in-plane) rotation, and a spring with 
rotational stiffness K w resists dw/dx (out-of-plane) rotation. 
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Assume that springs of the two types mentioned are present, with stiffnesses, placement, and gen- 
eralized deformation given by K:,y, «,;(«= 1, 2, ... ,N) for the first spring type, and 

K'w, y, w, r , (r — 1, 2, ... , R) for the second spring type. The elastic strain energy U associated with these 
springs us given by 


"-X 

n = 1 


~Y « - «>’) 2 l.-o i + X "T « - «>') 2 L 


r = 1 


OX 


(3.3.1) 


where it is assumed that the springs are unloaded at the rest state of the imperfect structure, 
{«} = {o°}. Using the expanded form for displacements along with the functional form of the various 
shape functions (see equations (2.4.1), (2.6.7), and (2.7.4)), equation (3.3.1) is evaluated and expressed 
as 


U = Constant + ( qq - q,q°)(F“ + F”) 


(3.3.2) 


where summation over the repeated indices is implied, the "Constant" term has no dependence on the 
modal amplitudes, and 


N 

[i + < -i a i,- 

n - 1 

F ti = C 1 + ( - 1 A ( X ) M>;) I y r 


(3.3.3) 


where functions £,{y) and <t>,{y) are the y-dependent portions of the eigenfunction components u,{x, y ) and 
W *( X ' y)’ respectively (see equation (2.6.7) ), and/' denotes differentiation of function /with respect to 


The elastic strain energy of the springs given in equation (3.3.2) is added to the total potential en- 
ergy of the structure. To obtain the stationary total potential energy condition governing equilibrium, the 
following derivative is required: 


dU_ 

a?, 




(3.3.4) 


Comparing equation (3.3.4) with equation (2.8.29), it can be seen that to account for the elastic support 
discussed here, each ij* coefficient C y and Cj of equation (2.8.29) is augmented by the quantity 

(f'l* + Fij). 


3.4 Constraint of the End Displacements 


r ^ le buckling analysis of VIP ASA simulates an infinite-length prismatic structure supported at 
uniform intervals along the length against transverse buckling displacements. It is the conventional as- 
sumption that the buckling eigensolutions will model closely the behavior of a finite-length structure if 
the halfwave lengths of the buckling modes are selected to be integer fractions of the length of the panel. 
The ends of the finite-length structure are assumed to be simply supported (see Section 2.3.2). In an 
effort to extend the method of NLP AN to structures with a greater variety of end-support conditions, a 
method has been developed for imposing constraints on certain generalized displacement components 
at the longitudinal ends of the structure. By using appropriately selected sets of buckling modes as shape 
functions and applying the generalized displacement constraints, a greater variety of end-support condi- 
tions can be simulated. 

Two different types of displacement constraints can be imposed at the ends of a panel in older to 
simulate three different types of support conditions. With the first type of constraint, the axial dis- 
placement u at specified points on the end of a panel is required to be equal to the effective axial dis- 
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placement of the panel end. This type of constraint is used to simulate either eccentric application of the 
end load (taken with respect to the neutral bending axis of the panel) or clamping of the panel end 
against wide-column buckling motion. With the second type of constraint, the slope dw/dx is constrained 
to be zero at a set of specified points at the ends of the panel. This type of constraint is used to simulate 
clamping of individual plate elements against out-of-plane rotation. The derivation of the constraint 
equations for both types of constraints is presented in this section. The incorporation of the constraints 
into the general procedure of obtaining equilibrium solutions is described in Section 3.5. 

The individual VIPASA buckling modes used as shape functions are not suitable for modelling the 
response of a structure having the modified forms of end support mentioned above. However, using a 
set of appropriately selected shape functions in conjunction with a set of generalized displacement con- 
straints, the desired types of response can be simulated. The point constraints are satisfied exactly, so 
in order to avoid having an over-constrained system, the number of point constraints must be fewer than 
the number of shape functions affected by the constraints. To simulate a stiffened structure with clamped 
ends, it is recommended to constrain the gross rotation of the ends using a couple of strategically placed 
axial constraints. The zero-slope constraint is suitable for simulating a simple rectangular plate with 
clamped ends. It is unwise to mix the two types of displacement constraints because the number of shape 
functions needed to provide meaningful results becomes excessively large. 

3.4. J Constraints on axial displacements. Locations at which the axial displacements are con- 
strained are specified in terms of points y k on the cross-section, where k is the index number of the point 
on the cross-section, and y is used here as a generalized transverse in-plane coordinate (ie. y signifies 
both a particular plate strip and a point along its local y-axis). The displacement u is constrained at these 
cross-sectional points, at both longitudinal ends, to be equal to the effective axial displacement of the 
respective panel end. In the current implementation, when axial displacement constraints are used it is 
required that displacement control of the in-plane loading be used (CONTRIV'D' in Table 1). With this 
limitation, the effective change in length Au, ff is given by 

Au tff =XAut (3.4.1) 

where X is the displacement-control parameter, and Auc is the value of Am corresponding to the linear 
response to the unit in-plane load system. 

An additional contribution must be added to the general expanded displacement form given in 
equation (2.4.1) or equation (2.8.1) for the following reason. The axial component of the primary shape 
functions (VIPASA eigenfunctions) denoted u,{x, y), has the following functional form on each plate strip 
of the structural model: 


t . . m.nx 

M, = %fy) cos -i- (1=1, 2 M) (3.4.2) 

where m, is the longitudinal halfwave number of the mode shape. The effective neutral bending axis 
of the panel is determined by the zeros of the functions £,< y ) (one function on each plate strip), since 
these zeros identify points on the panel end about which the end rotates during buckling. If the axial 
displacement of arbitrary points of the panel ends are to be constrained to zero in such a way as to 
simulate a different, eccentric, line of load application, then a corrective displacement contribution must 
be included which effectively shifts the zeros of the functions 2j,(y) from their original locations to new, 
specified, locations. For cases in which the modified function Ui = n + BiUe is used (see Section 2.6), the 
contribution B,Ub is included in evaluating the axial displacements at the ends. 

The corrective contribution to displacements {u} is taken to be q «{ml), where {iq.) is the linear 
unbuckled response of the panel to the unit in-plane load system (see Section 2.5), and where q. is an 
initially unknown amplitude parameter. The expanded form of the displacements is now given by 

(m) = X{u L } + q e {u L ) + ?,{«,} + qjqjlUij) (3.4.3) 

where summation over repeated indices is implied, and 
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(3.4.4) 


{«,} = {«,} + £,{%} 

Note that {u, y } which appears in equation (2.8.1) is replaced here by {«*}. This is because the restriction 
to displacement control (CONTRL='D' in Table 1) enforced here results in the equivalence of {«,,} and 
{ Uij) (see Table 4). In anticipation of later developments, parameter q, is represented as a series in the 
modal amplitudes, q,: 

4‘ = L fii (3.4.5) 

where summation over i is implied, and L (/= 1, 2 M) are constant coefficients. The displacement 

now are written as 


{“) = *•(«*,} + + qflj{Ujj} 

(3.4.6) 

where 


{“.) = («,} +B,{u b ) + L t [u L ] 

(3.4.7) 


The condition at each constrained point is stated as follows: the displacement u, minus the initial 
value associated with imperfections u°, minus a second value {^associated with a stress-producing forced 
rotation of the panel ends, is equal to the effective axial displacement of the panel end: 


(u k -(u 0 ) k -(i/) k ) = Xu L 


(£= 1,2 K) 

(x = 0, L) 


(3.4.8) 


where m* - u(y*), and K is the number of cross-sectional stations at which displacements are constrained. 
It is noted that because of the functional form of 
“vL=o = (3.4.9) 


The constraint equations (3.4.8) are evaluated in terms of the form for u of equation (3.4.6) and the 
simplification of equation (3.4.9) is applied, providing the following condition: 

(<7« -<&- + b Vb + LiU[) = 0 fin n'"’^ (3.4.1 0) 


where the parameters q{ are modal amplitudes used to specify forced rotation of the panel ends Using 
the functional form for u, of equation (3.4.2), 


I x = 0 — SiO) 

«il x=z , = (-ir^ l <y) 


(3.4.11) 


If m, is odd, then the corresponding values of B, and L, may be non-zero. However if m, is even then the 
corresponding value B, is zero (see Section 2.6), and it is stated in anticipation that the corresponding 
value L, is also zero. Equation (3.4.10) evaluated at x = 0 is subtracted from equation (3.4.10) evaluated 
at x = L, and the substitutions of equation (3.4.1 1) are used, to obtain 

X ( < ?<- < ?-°-<7f)(-2^ + B 1 Au fl + L 1 Au i ) = 0 (£= 1, 2 K) (3.4.12) 

‘ = 01.02,... 

where = ^{/), parameters o* ... are the index numbers of modes for which m, is odd, and A Uc and 
Aue are the contributions to A u (the change in panel length) due to Uc and respectively. 

The coefficients L, are determined by arbitrarily requiring that the first equation (k = 1) of equations 
(3.4.12) be satisfied identically for any value of the modal amplitudes q it The quantity inside the right 
set of parenthesis is set to zero for each value i % and parameters Li are determined to be 
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(3.4.13) 


L = B Au b (* = 1 ) 

‘ A u L ‘A u L O' = o lt o ^, ...) 

If AT is greater than one, then the remaining K — 1 non-trivial constraint equations are obtained by sub- 
stituting h of equation (3.4.13) into equation (3.4.12). The final form obtained is 

(<7i ~ <7 i - = 0 (* = 2, 3, ... , AO (3.4. 14) 

where 

r _ f ( “2£? + B, A« B + L, Au^) (r = Oj, ...) 

*■“ \o C *-«„<*.-) 0.4.15) 

where £\ f e 2 , ... are the index numbers of modes for which nu is even. 

Next, equation (3.4.10) evaluated at x = 0 is added to equation (3.4.10) evaluated at x = L t and the 
substitutions of equation (3.4.1 1) are used. Noting that B, and L, are both zero for i = e u e 2j the result 
is 


X ^-<7.°- ^(2^ = 0 (I: =1,2 AO (3.4.16) 

= ... 

These equations are distinct from equations (3.4.14), but can be expressed in the same form: 

(<7; ~ Qi ~ <Ii)Eia = 0 (k= 1,2 A0 (3.4. 17) 

where 


f 0' = e,,e 2 , ...) 

\ 0 0 = O], £>2, ...) 


(3.4.18) 


3.4.2 Constraints on out-of-plane rotation. The locations at which out-of-plane rotation is con- 
strained are specified in terms of points y* on the cross-section, where k is the index number of the point 
and y is used here as a generalized transverse in-plane coordinate (ie. y signifies both a particular plate 
strip and a point along its local y-axis). The rotation dw/dx (referred to a local plate-strip coordinate 
system) is constrained at these points on both longitudinal ends of the cross-section. At each end of the 
panel, K constraint conditions are imposed: 


w* = (w")* + (w,{)* 


(*=1,2 A0 

(x = 0, L) 


(3.4.19) 


where w , i is the slope at the cross-sectional station y 4 , w,° corresponds to the imperfection shape of the 
panel, and w/ z corresponds to a forced rotation of the ends. When equation (3.4.19) is evaluated in terms 
of the expanded form for w and the functional forms for w, and w y , the constraint equations at x = 0 have 
the form 


(<7. -q°- = 0 (*=1,2 A0 (3.4.20) 

and the constraint equations at x = L have the form 

Oft -Q° 1)'"’ = 0 (* = 1,2 A0 (3.4.2 1 ) 

where summation over i is implied, m, is the halfwave number of the i* buckling mode shape, <t>,(y) is 
the y-dependent function corresponding to w,(x, y) (see Section 2.6), <|>* = 4>i(y*), and { cf } is a set of model 
amplitudes used to specify a forced rotation of the panel ends. There are 2 K constraint equations; how- 
ever if halfwave numbers m, are either all even or all odd, then the two sets of constraint equations 
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(3.4.20) and (3.4.21) are equivalent, and there are only K independent constraint conditions. Each of the 
K or 2 K constraint equations can be written in the form 

(<7. ~<h ~ </)£*, (3.4.22) 


where 

E ki = rntfi (k = 1 , 2 , ... , K) 
and, if both odd and even values of nu are present, 

= (k = K+l,K + 2 2/Q 


(3.4.23) 

(3.4.24) 


3-5 Solution Procedure with Thermal Loading and Displacement Constraints 

This section describes modifications to the NLP AN problem formulation and solution procedures 
necessary to accommodate thermal loading and/or displacement constraints. The nature of the thermal 
loading, and the associated additions to the assumed form of the displacements, are discussed in Section 
3.2. The nature of the displacement constraints, the displacement-constraint equations, and the associ- 
ated modifications to the assumed form of the displacements, are all discussed in Section 3.4. In this 
section, the general form of the displacements which accounts for the new features is used to form a 
new total potential energy expression, and the constraint equations are incorporated by using the 
Lagrange multiplier method. The adaptation of the advanced solution strategies to algebraic equations 
containing Lagrange multipliers is also discussed. 

35.1 Energy functional with Lagrange multipliers. When both thermal loading and axial- 
displacement constraints are imposed, the displacements have the expanded form 

{«} =Mu l ] +Y{«/-} +qi{Ui) +q,<}j{Uij} ,i,j= 1,2,3,... (3.5.1) 

where summation over repeated indices is implied, y is the thermal-load control parameter, (u r } is the 
unit thermal response discussed in Section 3.2, and 

{«,-) = {«;) +B,{« b } +L i {u L ] 

{«(,) = {«,>} - {u L }) (3 ' 5 - 2 > 

where the form of [it,,] is explained in Section 2.7, coefficient B, (i= l, 2 , ... ) is defined in Section 2.6, 
and the presence of coefficient L, (i = 1 , 2, ... ) implies that constraint equations may accompany the 
problem statement. As stated in Section 3.4, axial displacement constraints are permitted only when 
displacement control is used (CONTRL='D' in Table 1), and this has the implication that 

L < = 0 0 = 1 . 2 , ... ) -or- A^s 0 (i',; = 1 , 2 ,...) ( 3 . 5 . 3 ) 

The expansion of the mid-surface mechanical strains given in equation (2.8.3) is modified here to 
include one new term and one altered term: 

{e m } =X{e i ,} +y{e”) + (q i -q‘){e i } +(.qflj-q°qj){i i j\ + ... ( 3 . 5 . 4 ) 

where y is the thermal-load control parameter, {e"} is defined in equation ( 3 . 2 . 4 ), and 

{£<} = (£,} +£,{%} +L i {t L ) ( 3 . 5 . 5 ) 

The expansion of the in-plane stress resultants now has the form 

{N} = \{N l ) + y [N t ] + (< 7 , - < 7 ,°){A/,} + (qflj - qfqJ){Ny\ + ... (3.5.6) 

where {/V r ) is defined in equation (3.2.7), and 
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{A/,} = {^}+B ( {/^J+L I {Ai} 


(3.5.7) 


(3.5.8) 


The total potential energy n of equation (2.8.22) is evaluated, using the expressions for (e") and 
{/V} of equations (3.5.4) and (3.5.6) in place of the original expressions given in equations (2.8.3) and 
(2.8.7), respectively. The resulting expression for 7 t is written as 

k = Constant 

+ «<xcf + TCf - pc? - ,/C? - q^cj, - qtfqfcju) 

+ - PC? + 1 C> - qtCfj - qtqfC?) 

+ W/lACqk + tcjjt + C?J - qfclt) 

+ + yclu + Cjl + -|- C~) + CKqf) 

where the "Constant" terms are those which do not depend on the modal amplitudes. The coefficients 
appearing in the above equation follow the definitions given in equations (2.8.25) except as noted here 
Four coefficients have additional contributions, given by 

ACf = L,c£ 

AC? = L,cf + Lft + L,Lft + (LjBj + B^)C L b 
^ = L k cfj 
A c£* = Z*c£* 

and two coefficients, not present in equation (2.8.24), are given by 

cJ = c] + B t C T B +L t C T L 
Cjj = cjj + Ay<cJ - C[) 


(3.5.9) 


( 3 . 5 . 10 ) 


Seven primitive coefficients (Cl, Cb, Cj, Cj , Cj Jt Cj A , Cjju ) appear in the above equations, and these are 
defined as illustrated here by two examples: 


P 

c ^X(J/ { ^ }7 ’ {e v ))dL4 ) 


(3.5.11) 


p= i ' ' p 

where the notation used is the same as in Section 2.8.3. 


The constraint equations for the different types of generalized displacement constraints all have the 
same form, visible in equations (3.4.14), (3.4.17), (3.4.22). It is assumed that there are a total of N 
constraint equations, written collectively as 

(< 7 > ~ <? ij ~ dj) E nj = 0 (« = 1 , 2 N) ( 3 . 5 . 12 ) 

where summation over j is implied. The constraint equations are incorporated into the problem statement 
using the Lagrange multiplier method. A new functional n is formed: 

n = 7C + r„(< 7 y - q° - q{)E n j (3.5.13) 

where T. (n = 1,2, ... . AT) are die Lagrange multipliers, and where summation over repeated indices is 
implied. The equilibrium condition for the constrained system is obtained by setting the first variation 
of the new functional to zero: 
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(3.5.14) 


5n = 


an 


a?, 


&7,+ 


an 

ar. 


= o 


sr„ 


Because each value 5q, or 8r. is both independent and arbitrary, equation (3.5.14) gives rise to the fol- 
lowing (M + N) independent conditions: 


an 

a?, 


+ q/Xcjj + ycjj + pcf + Cy - q° k C*j - q° k q?c“) + r„F m 

+qM^+ycl>+c, Jk -q?c! Jk ) 


+ qflkQ&c'iju + y cfjki + Cijki) 
= 0 (i= 1, 2, ... ,M) 


an 

dr„ 


= (qj-qf-q!)E nj 
= 0 (n= 1,2, 


AO 


(3.5.15) 


(3.5.16) 


where the coefficients of equation (3.5.15) are all defined in equation (2.8.30) except for the following: 


Fin = Eni 

cj=cj 

At rr 

rf = 2C • 

At 'r t 

Cqk = cjj k + 2CJ& 

cjw = 2(^ + 5,) 


(3.5.17) 


Equations (3.5.15) and (3.5.16) govern the equilibrium solutions which are sought. There are now 
three generalized load-control parameters, X, y, and p. The modal amplitudes q, and Lagrange multipliers 
r„ constitute (Af + AT) variable parameters to be determined in obtaining each equilibrium solution. 

3.5.2 Modifications to the solution strategies. The solution procedures described in Section 2.9 and 
Appendix D require some additions and modifications to accommodate the additions to, and modifica- 
tions of, the nonlinear algebraic equations discussed in the preceding section. First, the new parameters 
appearing in equations (3.5.15-16) must be normalized. 

The thermal load parameter y is normalized by a reference value y^ so that the normalized load 
parameter y is given by 

y-yhrtf (3.5.18) 

The value y, €/ is selected to be the maximum-amplitude target value for y used in specifying the load 
range or ranges over which a nonlinear analysis is to be performed. 

The Lagrange multipliers T. are normalized by reference values R, so that the normalized Lagrange 
multipliers T. are given by 

r„ = r n /R n (/i =1,2 AO (3.5.19) 

The values R* are selected as follows. For each value of n the value of i is identified which corresponds 
to the largest-amplitude coefficient F„ for all t, where appears in equation (3.5.15). For the selected 
value of i, the following expression is set to zero: 

Cu + r^,„ = 0 (3.5.20) 
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where no summation is implied. Equation (3.5.19) is solved for T. and that expression is used in equation 
(3.5.20) to get 

Q + W^ = 0 

Parameter R , is determined by setting r„ to unity in equation (3.5.21), providing 
R n = -CJF I „ 

The normalized parameters y and T, are incorporated using the same approach as was described in 
Section 2.9.1 for incorporating the normalized parameters X and 0. The nonlinear algebraic equations 
(3.5.15) are expressed in terms of the normalized parameters, and the resulting equations are identical 
in form to equations (3.5.15), but with normalized parameters y and T. replacing y and T., respectively, 
and with the constant coefficients having been transformed in the obvious manner to account for the 
definition of y and T.. Without assigning new notation, it is assumed from here on that the parameters 
and coefficients which appear are the normalized ones. 

Following the general procedures described in Sections 2.9.2 and 2.9.3, the form of the nonlinear 
algebraic equations is simplified after specification of the modal imperfection amplitudes {q°} and col- 
lapsing of the load parameters X, y, and 0 into a single load parameter A. Equations (3.5.15) and (3.5.16) 
are converted to obtain equations of the following form: 

(C,- + ACf) + q/Cij + A C*j) + qj q k (C ijk + ACj*) 

+ + A cjj U ) + r n F in = 0 (i- 1, 2, ... , M) (3 - 5 ' 23) 

E n + Qjpnj = 0 (n = 1, 2 N) (3.5.24) 

where the coefficients of equation (3.5.23) are analogous to those of equation (2.9.12) (with the tilde 
dropped from the notation), and where 

E n = ~ (q° + <Jj)E nj (3.5.25) 

where summation over j is implied, and it is assumed that the modal amplitudes {cf} which determine 
the forced end rotation have been specified. 


(3.5.21) 

(3.5.22) 


Equations (3.5.23) differ in form from equations (2.9.12) by the presence of the terms with T., and 
by the accompanying constraint conditions of equations (3.5.24). Equations (3.5.23) and (3.5.24) can 
be converted into a set of M+N equations with the same form as equations (2.9.12): 

(Cj + AC, ) + qfCij + AC”) + qjq k (C i]k + AC”*) 

+ qm&iju + AC ijki ) = 0 (|'=1,2 M + N) (3.5.26) 


where repeated indices are summed over the range 1 through M + N , and 

Qi = E n • i=M + n («= 1,2, ... ,A) (3.5.27) 

The only non-zero coefficients in equations (3.5.26) with subscript values i, j, k, or /, greater than M , are 
given by 


c ij = F m ’ j = M+n 


0 =1,2, ... ,M ) 

(n= 1, 2, ... ,N) 


E i~ E n > i = M + n (n= \,2, , N) 


(3.5.28a) 

(3.5.286) 


Cij = E nj . i = M+n 


(*=1,2 N) 

0= 1.2 M) 


(3.5.28c) 
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Many aspects of the solution strategies discussed in Section 2.9.4 and Appendix D can be applied 
directly to equations (3.5.26). However while there are M + N variable parameters (excluding the load 
parameter) in the M + N equations, there are only M-N independent geometric variables (modal am- 
plitudes) because of the N constraint conditions of equations (3.5.24). Thus, the eigenvalue problems 
which are solved as a part of the solution strategies described in Appendix D have only M-N 
eigensolutions, and so some modifications must be made to the solution strategies. 

Appendix D Section D.3 concerns the control of solution branching using Thurston's method. 
Equations which appear in that section are modified here to accommodate the use of Lagrange multi- 
pliers. The M + N equations represented by equations (3.5.23) and (3.5.24) are written symbolically as 

f,{q, T, X) = 0 (i = 1, 2, ... , M + N) (3.5.29) 

where the over-bars signify vectors. Let (g,r , X) be a starting solution (a known exact or approximate 
solution to equation (3.5.29)) and let (q + £, r + £, A, + 5) be a solution which is sought near the starting 
solution, where the increments to the parameters are small compared to the parameter values The new 
solution satisfies 


^<<7 + tr + C,X + 5) = 0 (i= 1,2,... ,M + N) (3.5.30) 

The first M of equations (3.5.30) are expressed in the expanded notation of equation (3.5.23), and terms 
are regrouped based on their order in the incremental parameters: 

(D, + 5 D, 6 ) + 5/Di, + soj) + (yU(D„, + SD?J + %£UD, tU + ic^) 

+ t/i = 0 (,= 1.2 M) <3 - 5 - 3 ') 


where 


D,=m r,X) 


(3.5.32) 


and the other coefficients are defined in Appendix D. The last N of equations (3.5.30) are expressed in 
the notation of equation (3.5.24): 

( £ n + + ^„, = 0 (/i = 1 , 2 N) (3.5.33) 

Equations (3.5.24) are linear in the modal amplitudes, and hence these equations are always satisfied 
exactly, even in an iterative solution procedure. Thus the quantity in the parentheses in equation (3.5.33) 
is known to be zero, and the equation reduces to 


^ = 0 (*=1,2 N) 


(3.5.34) 


The eigenvalue problem for the constrained system here, corresponding to Appendix D equation 
(D27), is 




(*= 1 , 2 , ...) 


(3.5.35) 


where 8* are eigenvalues, and [>|T*] r are eigenvectors with components coircsponding to the parame- 
ters UlU • A senes of operations is used in NLPAN to reduce equations (3.5.35) to a system of 
M-N equations which can be solved to obtain M-N eigenvalues and eigenvectors. The details of this 
process are tedious, but straightforward. The first step is to reorder the components of (0 1 } so that the 
first N columns of the reordered matrix [£] and the first N rows of the reordered matrix [£] form 
non-singular square matrices. The last N rows of equation (3.5.35) are then used to express N compo- 
nents of {9 1 } in terms of the remaining M-N components. Similarly, the N components of {t*} can be 
expressed in terms of the M-N remaining components of {0*} and the eigenvalue 5*. A condensed 
eigenvalue problem of dimension M - N is thus obtained, and this is solved using conventional methods. 
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In Appendix D equation (D29), the incremental solution vector {^} is expressed as a series in the 
eigenvectors {0*}. The analogous equation for the modified case here is given by 

(35 36) 

The other aspects of the solution strategy described in Appendix D Section D.3 are not affected by the 
presence of Lagrange multipliers, provided that: i) the M + N equations of the form of equation (3.5 26) 
are used as the equilibrium equations (the Lagrange multipliers r, are treated the same as modal am- 
plitudes), and n) it is realized that there are only M-N eigensolutions. For example. Appendix D 

equation (D25) would be of dimension M + N whereas Appendix D equation (D31) would be of di- 
mension M-N. 

Appendix D Section D.2 concerns the application of the arc-length-control solution strategy. Once 
/•f so * ut ' on strate gy can be used for the most part without modification, assuming that equation 

(3.5.26) serves as the starting set of equations. However special treatment is required is in solving the 

eigenvalue problem of Appendix D equation (D20). The corresponding equation for cases with Lagrange 
multipliers is 66 

{ -^r } = {■§■} <*=1.2.») (3.5.37) 

where to* are eigenvalues, [<J>* 1 are eigenvectors corresponding to the variables [t|CY, and [/] is 
the identity matrix. Equation (3.5.37) has the same general form as equation (3.5.35), and thus the 
method used to reduce equation (3.5.35) to a problem of dimension M-N is used again to reduce 
equation (3.5.37) to a problem of dimension M-N. 
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4. IMPROVEMENTS IN THE COMPUTATION OF THE SECOND-ORDER DISPLACEMENT 

FIELDS 

In this section, consideration is given to alternate methods for computing the second-order dis- 
placement fields. This is done in an attempt to rectify some dilemmas and discrepancies which arise 
when the method described in Section 2.7 is used without modification. First, the troublesome aspects 
of the second-order fields are identified. Next, the beneficial effects of requiring the second-order fields 
to be orthogonal to the buckling mode shapes are discussed. Several possibilities for new or modified 
approaches to computing the second-order displacement fields are outlined. Finally, the approach used 
in the current implementation of NLP AN is described. 

4.1. Problems with the Second-Order Displacement Fields 

4.4.1 Boundary conditions at the longitudinal ends. The buckling modes {«,} are the primary shape 
functions in the NLP AN analysis, and their harmonic dependence on x determines the effective boundary 
support condition at the longitudinal ends of the structure. The functional form of {«,} of equation (2.6.7) 
suggests that boundary support is in place which guarantees that 

v = w = 0 at x = 0 , L ( 4 . 1 . 1 ) 

The linear prebuckling solution ( 14 .} does not rigorously satisfy equations (4.1.1), because of in-plane 
Poisson expansion, and the coupling of in-plane and out-of-plane displacements between adjoining plate 
strips. In addition, and of primary concern in the current discussion, the two contributions to both v 0 and 
wu of the second-order fields have a cosine dependence on x (see equation ( 2 . 7 . 4 )), and thus it is not 
guaranteed that v, y and w v are zero at x = 0, L. Indeed, it was observed in Ref [7] that if left unchecked, 
some contributions to selected functions {«,>} grossly violate the boundary conditions at the longitudinal 
ends. 

The functional form of { 14 } was selected to accommodate the non-homogeneous terms in the gov- 
erning differential equations (equations (2.7.1) ) so as to permit the equations to be expressed in terms 
of separated variables. Boundary conditions at x = 0, L were not considered in selecting the functional 
form of Uij, but justification for the functional form is found by considering the behavior of simple rec- 
tangular plates. 

In the postbuckling analysis of simple rectangular plates using von Karman plate theory, where 
w(x, y ) is represented, as it is here, as a series of terms of the form <j>,<y) sin (nunx/L), the functional forms 
for Uifcx, y) and v,/x, y) of equation (2.7.4) are the appropriate forms to allow the in-plane equilibrium 
equations to be satisfied exactly for any arbitrary set of terms in the series for w. (The in-plane boundary 
conditions at x = 0, L corresponding to % and v,y are u„ = 0 and AX, = 0.) The in-plane displacement terms 
un and v,y are essential for obtaining accurate solutions for a plate undergoing significant postbuckling 
deflections. 

For a linked-plate configuration, the coupling of the transverse displacement components v and w 
between adjoining, non-coplanar plate strips suggests that component should have the same harmonic 
dependence on x as component v ;; . It is concluded that with the use here of global functions, the boundary 
condition of equation (4.1.1) is too restrictive. It seems more appropriate to require that at the longi- 
tudinal ends of each plate strip the out-of-plane displacements must be small, and are permissible only 
to the extent that in-plane displacements must be accommodated and displacement compatibility at the 
node-lines must be enforced. The boundary value problem described in Section (2.7) requires inter- 
vention if the qualitative boundary condition just described is to be enforced. 

4.1.2. Load-dependence of the second-order fields. The differential equations (2.7.1) which govern 
{%} were are obtained from an expansion of the plate equilibrium equations in terms of the modal 
amplitudes. The appearance of the load parameter X in equation ( 2 . 7 . 1 ) suggests that the functions {14 } 
are load-dependent; this is an undesirable quality from the standpoint of computational economy. A 
second complication related to the appearance of X in the differential equations is the existence of sin- 
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gular values of X for which the amplitude of (u,,) becomes unbounded. Understanding the reason for 
these singularities gives insight as to how {«,y} might be computed differently, therefore the cause of the 
singularities is discussed here. 


When the buckling equations (equations (2.6.4), (2.6. Ib-c) ) are expressed in terms of separated 
variables, the following homogeneous ordinary differential equations governing the function {£•} are 
obtained: 


{ W+CA+Ci n/ 1 r 0 
D&' + D* ryr + Dtfi l + X*I 

E^,"" + E&" + £3<t>, J [ far + 



(4.1.2) 


where primes denote differentiation with respect to y, and the sub- and super-scripted coefficients C, 
D, and E depend on the halfwave number nu. Similarly, when equations (2.7.1) are expressed in terms 
of separated variables, the following nonhomogeneous ordinary differential equations governing the 
functions {^} are obtained: 


{ -C&J-C&t + Ci tW 1 f 0 If 1 

-D&v' + D^’ + D^j V + W D 3 V> f = 1 G(U,M$,}) > (4.1.3) 

E^'"' + E2^/' + J t + J f //({£,]. {£,}) J 

where the sub- and super-scripted coefficients C, D, and E are the same as those in equations (4.1.2), 
except that they depend the halfwave number m of equation (2.7.5) instead of nu. It can be seen by 
comparing equations (4.1.2) and (4.1.3) that if X in equation (4.1.3) is selected to be an eigenvalue for 
a buckling mode with the halfwave number m, then the left-hand side of equation (4.1.3) (when ex- 
pressed in terms of unknown coefficients for the functions {£«->}) will be singular. (This fact was pointed 
out to the author by Prof. S. Sridharan of Washington University at St. Louis.) Thus for some ranges 
of values of X, the functions {«,>} take on large amplitudes and have the approximate shape of buckling 
modes, except that the phase of v,y and w :j causes the maximum displacement amplitudes to be at 
x = 0, L, resulting in a gross violation of boundary condition. In these situations the second-order 
functions {«,>} are reflecting instability-related response which is already represented by first-order 
functions {«,}. 


4.2. Implications of Imposing Orthogonality Between [u^} and 

In the literature of the classical perturbation approach to the analysis of structural stability, it is 
stated that second- (and higher-) order displacement fields should be orthogonal to the buckling modes 
(Ref. [22]), although it is not necessarily clear how this orthogonality is to be enforced. In this section, 
the implications of requiring orthogonality between the functions {«,} and {«,>} are studied. Evidence is 

offered that the problems discussed in the preceding section are alleviated by enforcement of this con- 
dition. 


Abbreviated notation is introduced here. Define operators L and N such that L(N) and N(N,u) are 
given by 


L(A0 = 


l f 

w J = { 

J l 

r Ni (N y , U ) i r 

= <{ N 2 (N x ,v) j. = <{ 

l N 3 (N,w) J [ 


N +N 

1 T z»x ' 1 y xy*y 
^y>y 

^X*XX ^^Xy'Xy ^y*yy 


} 


(A fyU,y),y 

(N x v, x ), x 


(N X W, X + A^W, y ), x + (A^W, X + NyW,y),y 

Define the inner products <u, L(/V)> and <u, N(N,u)> by 


} 


(4.2.1) 

(4.2.2) 
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p 

L(AO > = ^ I J C" ^i(AO + vL^AO + w L^{M)~\dA J 
N(N,u)> = Su \_uN x (N y ,u) + vN 2 (N x ,v) + wN 3 (N,w)ldA J 


<«, N(N,u)> = luN { (N y ,u) + vN 2 (N x ,v) + wNjiN.wfldA J 

Using the abbreviated notation, the plate equilibrium equations are given by 

L(N) + N(N,u) = 0 (424) 

the buckling equations are given by 

L(N,) + X,N(N l ,u) = 0 (4.2.5) 

and the equations governing {«,>} (equations 2 . 7 . 1 ) are written as 

L W/> + y “j) + y W(ty «i) + XN(N l , UiJ ) = 0 ( 4 . 2 . 6 ) 

where it is noted that in equations (2.6.4) and (2.7.1), it is assumed that Ni(Nyt ,, u)= 0 . 

The orthonormality relationship satisfied by the buckling modes is derived in Appendix B and given 
in equation (BIO). It is written here as 


<u r N(N l , u,)> = - h tj a t 


(4.2.7) 


where 5,y is the Kroniker delta functions, and a, is an arbitrary normalizing constant. The condition ex- 
pressing orthogonality between the second-order displacement fields and the buckling modes is given 
in Appendix B equations (B 12) and (B 14) in two alternate forms: 


<«*. N(N l , u tJ )> = 0 
<u ijt N(N l , « t )> = 0 


(4.2.8a) 


Enforcement of the above orthogonality conditions would cause the second-order fields to have 
components w tJ (in the local reference systems) which are closely in accordance with the goal of having 
minimal transverse displacements at the panel ends. This is because the large transverse displacements 
which characterize the buckling modes would be suppressed in the second-order fields by the 
orthogonality condition. The second-order fields would thus be truly second-order in character, absent 
of the displacement contributions already available in the family of buckling modes. 

It is also contended here that enforcement of the orthogonality condition would cause the second- 

OI ? e u fi ! ,dS 10 56 inde P endent of ^ load parameter. The term X, N(N L , u) of equation (4.2.5) is the term 
which dnves the buckling instabilities associated with the mode shape {«,}. Function {«,,} is now 
orthogonal to all the buckling modes and thus should not represent an instability-driven displacement 
field; therefore it would seem reasonable that the term XN(N l , u,,) of equation (4.2.6) should have little 
effect on the solution for {«, y ). 


A more theoretical basis for this argument is established by considering terms in the stationary total 
potential energy expression. The total potential energy of the structure is expressed in equation (2.8.24) 
as a power series in terms of the modal amplitudes, which can be expressed symbolically as 

It = JCo + 7C] + ^2 + 7^3 + H 4 + ... (4 2 9 ) 

The equilibrium condition of equation (2.8.28) can similarly be written as 
811 = Stcj + 5 ^ + 5tc 3 + 5 jt 4 + ... 

= 0 (4.2.10) 
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It is noted that while the boundary-value problems governing the functions {«(.}, {«,}, and {u, y } were 
determined by expanding the plate equilibrium equations and the associated boundary conditions, iden- 
tical boundary-value problem statements for the first two functions can be obtained by evaluating 
5jii = 0 and 87 c 2 = 0, respectively. The relationship between 8k 3 and the boundary-value problem used in 
Section 2.7 to obtain {«, y } is more complex. 

Appendix C contains an evaluation of 8 tc 3 in which the following form is obtained (equation (C7)): 

6 jc 3 = ( - <«,, lUN Jk ) + -1 N(N k , uj) + \ N(Nj, uj) + XN(N l , «,*)]> 

' z z (4.2.11) 

-2 (k-\ k )<u i j,N(N L ,uj)>} 


In this equation, the weighted term inside the first inner product is the term set to zero in the differential 
equations (4.2.6), and the second inner product is the quantity which is set to zero in the orthogonality 
condition of equation (4.2.8b). If equations (4.2.8a-b) are enforced, equation (11) becomes 

8% = & 7 ,< 7 y< 7 * ( - <u„ lUN jk ) + -i- N(N k , uj) + y N(Nj, u t )]> J (4.2. 1 2) 

The load-dependent term XN(N l , uj) has been eliminated in the expression for 8713 . This is offered as 
further evidence that the second-order displacement fields are load-independent when orthogonality is 
enforced between the functions {u, y } and {«*) . 


4J Computation of {«<,} with Orthogonality Imposed 

With the load-dependent term omitted, the differential equations (4.2.6) are written as 

Wij) + y A Wj, Mi ) + N(N h uj) = 0 (4.3.1) 

A consequence of enforcing the orthogonality condition of equation (4.2.8) is that equation (4.3.1) can 
not, in general, be solved exactly over the domain of the structure. It seems appropriate to devise a 
method which minimizes the error of equation (4.3.1). Three possible approaches to computing {«*} are 
discussed in the following, then the approach currently used in NLPAN is described. 


4.3.1 Least-squares approach. One solution approach is to use the Lagrange multiplier method to 
satisfy the orthogonality-constraint equations while minimizing the error of the field equations in the 
least-squares sense. Equation (4.3.1) represents three equations. Define three residual error functions 
associated with the ordinary differential equations governing {^J (equations ( 4 . 1 . 3 )), where the load- 
dependent term is now omitted: 

R u (y) = - C£" - C£ + C 3 r\' - F(y ) 

R v (y) = -D£ + D 2 t\" + D 3 r]-G(y) (4.3.2) 

R w (y) = £,<)>"" + E^Sf" + Eji j) - H(y) 


where subscripts are dropped here, and throughout the remainder of this section, from 2 ^, q my , and <|> a 
The constraint equation (4.2.8b), is evaluated to get the equation 


a= 1 


1 ' f 


, f nnc , m 7C 

J 0 (_ z; 


mux v . f .... n 

) sin( — - — ) > =0 


} 


(4.3.3) 


If (m, m») are (even, even) or (odd, odd) then the x-integral of equation (4.3.3) is zero. Otherwise, the 
y-dependent portion of the equation is satisfied independently for each value of a; this is written as 
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n + ^(y) <t>] dy j =0 

' p 


where 


(4.3.4) 


H*O0 = [ - A^( ~ )nj . 0*(y) = l- Nxt {^- )\ + A^*"] 


(4.3.5) 


by 


The second-order displacement fields are now determined by minimizing the functional W, given 


'■|(f 


l(R u f + (R v ) 2 + (R w ) 2 ]dy 


); 4 (t 


C H *0) + 0*0) <t>] dy 


). 


(4.3.6) 


where T* (k = 1 , 2, ... ) are Lagrange multipliers, and where boundary conditions at the node-lines must 
also be met. "Hie functional W is minimized by setting to zero the first variation of W with respect to 
the functions q, r|, and <|>, and with respect to the Lagrange multipliers r*. 


Preliminary work has been done on developing a procedure for solving equation (4.3.6) using a 
finite-difference representation of the functions 2^, t), and <|>, such as is used in the existing solution 
method, described in [7]. In the solution procedure of [7], the generalized force-resultants and gener- 
alized displacements at the node-lines are isolated algebraically, so that the (homogeneous) node-line 
boundary conditions can be applied directly to the system of equations governing the finite-difference 
solution. It appears that this is also possible in a finite-difference solution of equation (4.3.6), except that 
the moment resultants at the node-lines do not appear in the resulting system of equations. It thus ap- 
pears that for each node-line which is unrestrained with respect to rotation, the zero-moment boundary 
condition must be enforced using a constraint equation which is incorporated into the functional W using 
an additional Lagrange multiplier. 


4.3.2 Subtraction of buckling mode shape contribution. In this approach, {«,,} is initially computed 
using the method described in Section 2.7, and is then modified by subtracting contributions in the 
shapes of buckling modes. The contributions are identified using the orthonormality condition for the 
buckling modes, equation (4.2.7). This equation is evaluated in terms of separated variables, and is then 
written in the following way, where new operator and inner-product definitions are introduced: 

p 

<{£,}. £({£,})> =X [ -( J r L + ww + », L (wr<iy] 3 

= -6 ipi 


where b { is a constant coefficient Similarly, orthogonality between {u,,} and fu*}, equation (4.2.8a), can 
be expressed in terms of the contributions for each value of a as 


<{£*), Z,({^ ay })> = 


= S c - ( - 5 r ) W J <i* n«, + ♦«,*'>' +«J ♦» 

p = 1 x y v o •'o 


p 

= 0 


(4.3.8) 


where {£<*>} is the modified function which satisfies the orthogonality constraints. 


pie unmodified function is expressed as the sum of the modified function and a series of 
buckling-mode contributions {£*}: 


k 


(4.3.9) 
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where coefficients c k are initially unknown, and 



The buckling modes {Jj*} (r= 1,2, ... ) used in equation (4.3.9) should, ideally, be limited to those with 
associated longitudinal halfwave numbers which equal the halfwave number |m| corresponding to 
the function {iUt,}. There are two reasons for this. First, displacement fields with different halfwave 
numbers are orthogonal by virtue of their longitudinal waveforms. Second, the relationship between 
transverse functions shown in equation (4.3.9) implies a relationship between full-field functions, and 
this relationship makes sense only for the case where equation (4.3.9) is applied on the basis of common 
longitudinal halfwave numbers. 

From equations (4. 3.7-9) it can be established that 

<{£*}. £({5a*})> = -cA (4.3.11) 

Coefficients c* (k= 1,2, ... ) are computed from the above equation, and then the function {£^} can be 
computed by applying equation (4.3.9). 

4.3.3 Direct suppression of displacements. In this approach, the method of Section 2.7 is used to 
compute {«,;}, but for selected functions {£,>} displacement constraints are imposed directly on the model 
during the finite-difference analysis. The placement of constraints is done on an intuitive basis; the 
general approach is to place constraints along node lines in a way which suppresses large transverse 
displacements, while still allowing in-plane expansion/contraction of plate strips. For example consider 
the blade-stiffened panel of Figure 3. It would be appropriate to enforce w = 0 (global) at node-line 
numbers 1, 3, and 4, and impose v = 0 (global) at node-line number 2. 

Based on results presented in [7], the suppression of displacements using this method is justified 
for fields {1^} with an associated halfwave number m of zero. It has also been found that suppression 
of displacements for m = 2 gives improved agreement of analytical results for postbuckling in column- 
like modes with predictions based on column theory. For large values of m, matching the boundary 
conditions at x = 0, L is less important than predicting the proper behavior away from the ends, so direct 
suppression of displacements seems inappropriate in this case. Base on results presented in [7], the vi- 
olation of boundary conditions at x = 0, L is less of a problem with the fields having large values m. 

4.3.4 Current approach used in NLP AN. In the current implementation, NLP AN uses a combination 
of the methods described in Sections 4.3.2 and 4.3.3. For long wavelength contributions to {u,y}, 

I m | <4, direct suppression of transverse displacements, as described in Section 4.3.3, is used. An au- 
tomated procedure positions the displacement constraints at selected node lines, assuming that the con- 
figuration is a conventional stiffened panel configuration such as a blade-, T-, hat-, z-, etc. stiffened 
panel. For unconventional configurations such as complex column sections, the automated procedure 
may place constraints inappropriately, so the specification of constraints needs to be done on a case- 
by-case basis. 

For values m = m, ± where m, is large and m, is small, or vice versa, the method of Section 4.3.2 
(subject to certain modifications) is applied. These fields, referred to in the literature as "mixed 
second-order displacement fields," are known to have a tendency to duplicate the shapes of buckling 
modes. Let m, be the (small) halfwave number for the global-buckling mode, and let mi be the (large) 
halfwave number for the local-buckling mode. For large m, and small m,. 

|m| = \m,±m g \ (4 3 12) 

=m, 
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In Section 4.3.2, it was stated that, ideally, for contributions {£,} to be subtracted from the field 

(£<.>}. In the NLP AN code, one or more modes {£,.} are already in use which satisfy rru^m,. These 
modes are used with the premise that they satisfy the approximate relationship | m | ; however this 

may or may not be the case in general, because modes are classified using the conditions m, < 3 and 
nu'tA. For the treatment of many stiffened panel problems, displacements u, and v r are zero for the fields 
of interest, so that matching of m, and m is not important. It is not known whether a mis-match in these 
values is detrimental to the accuracy of results for more unusual configurations. A more theoretically 
pure approach would be to carry along buckling modes which are not necessarily used as shape functions 
in the NLP AN analysis, but which match the halfwave numbers |m | encountered in computing the fields 
{«,,}, and which, thus, can be used to identify contributions to be subtracted from the second-order dis- 
placement fields. This latter approach has not been implemented. 
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5. MISCELLANEOUS CONSIDERATIONS 

This section includes brief discussions of some miscellaneous considerations regarding the use of 
the NLP AN analysis method and computer code. First, considerations in the design of geometric rep- 
resentations of stiffened panels are discussed. Next, the convergence of analytical results with respect 
to the finite-difference discretization and mode-set selection is discussed, and mode selection strategies 
are suggested for some common cases. Finally, the factors affecting computation time and computer 
memory requirements are discussed. 

5.1 Geometric Representation of Stiffened Panels 

For stiffened panels with multiple evenly spaced stiffeners, experience suggests that a unit-cell 
representation in NLPAN is generally preferable to a full, multiple-stiffener representation. (An example 
of a unit-cell representation is shown in Figure 5(c).) One reason is that with a unit cell, only a very few 
local-buckling modes (one to three) need to be incorporated in the analysis to allow the panel to take 
on various types of local deformation (for example, stiffener-web buckling, skin buckling, and flange 
buckling), whereas for a full-panel model, a much larger number of modes is required to accomplish the 
same thing. This is because the various stiffeners and skin bays in a multiple-stiffener model tend to 
participate to different degrees, and in different manners, in any given buckling mode. Another reason 
for using a unit-cell representation is that execution time and computer memory requirements increase 
with the complexity of the cross section. As a consequence, the number of buckling modes which can 
be incorporated as shape functions in an analysis decreases with the complexity of the cross section. For 
purposes of comparing the results of a unit-cell analysis with the results of a full-model analysis or test, 
it is suggested that the reference load values used in normalizing the various result sets be selected on 
the basis of a common axial strain value (assuming that the loading is uniaxial). 

5J2 Convergence Considerations 

5.2.1 Discretization of the cross section. The cross section of a configuration is discretized in order 
to perform both the finite-difference analysis described in Appendix C of [7], and the numerical inte- 
gration of the coefficient expressions such as those found in equations (2.8.26-27). Specifically, the 
y-dependent variables on each plate strip are evaluated only at a set of discrete, uniformly spaced points 
along the local y-axis. Increasing the fineness of the discretization improves the accuracy of results, but 
also increases computer memory requirements and increases program execution time. 

How rapidly the results converge with increasing fineness of the discretization depends on both the 
node-line boundary conditions, and the in-plane load conditions. For example, consider a square, simply 
supported plate subjected to a uniaxial compressive load N M . For the in-plane boundary condition given 
by BCVEC(2,IB)=2 in Table 3, the load N, is uniformly zero along the y-normal edges, whereas for 
BCVEC(2,IB)=3, the y-normal edge remains straight, and the average value N y is zero. Despite the 
seemingly small difference in these two sets of boundary conditions, in order to obtain similar accuracy 
in the predictions of postbuckling response for the two cases, the latter boundary condition requires the 
use of only about one third the number of discretization intervals as the former [7], For the former case, 
a minimum of 30 intervals is recommended, whereas for the latter case, ten intervals provides compa- 
rable accuracy. If the load axes are reversed for this problem so that uniaxial N, loading is applied, the 
number of discretization intervals required for a given level of accuracy is significantly greater than for 
either of the two cases just described. 

The minimum number of intervals allowed on any single plate strip is four. When using a unit- 
stiffener-cell representation of a uniaxially loaded stiffened panel, where symmetry conditions are im- 
posed on the skin at the edges of the cell, the use of a minimum of twelve discretization intervals for 
the skin to either side of the stiffener is recommended for local/global mode interaction problems. A 
convergence study is recommended as the best way to assure that a model is adequately discretized. 
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5.2.2 Mode set. The selection of VIPASA buckling modes used as shape functions in an NLPAN 
analysis affects the accuracy of the analytical results in two ways. First, the true qualitative response 
of the physical structure can be predicted analytically only if a suitable mode set is incorporated. Sec- 
ond, assuming that the appropriate family of buckling modes has been identified, increasing the number 
of modes used in the nonlinear analysis enables the computation of accurate results deeper into the 
nonlinear regime. 

In some situations, the set of suitable buckling modes can be selected using intuition. For example, 
the mode set for an axially compressed square, simply-supported plate with an edge-length a is given 
by sin(m7tx/a) x sin(mcy/a) (m, n= 1, 3, 5, ... ) where modes are added starting with the lowest values for 
m and n. For general NLPAN configurations, it is much less obvious what comprises an appropriate 
mode set, and unfortunately, the use of an inadequate set can cause errors ranging from erroneous 
stresses and strains, to the complete failure to predict some modes of response. 

For flat, stiffened or unstiffened panels subjected to in-plane loading, some guidelines for selecting 
modes are provided here. The guidelines are based on experience in modelling panels in which the 
buckling modes are classifiable as "global" or "local”. A global mode is characterized by a long wave- 
length and minimal distortion of the cross-section. A local mode is characterized by short wavelength 
buckling of plate strips in the structure, with significant distortion of the cross-section. Separate guide- 
lines are offered for symmetric structural sections and unsymmetric structural sections, because the for- 
mer can generally be modelled with fewer modes. ("Symmetric" refers to the initial geometry, not to 
the response.) 

The mode selection guidelines are presented in Table 5. In the table, the label (m,i) is used to 
designate the i* buckling mode in the infinite sequence of modes having the longitudinal halfwave 
number m, where the modes are ordered based on their eigenvalues. Label m is used to designate the 
longitudinal halfwave-number for the critical local-buckling mode. A modal-interaction analysis is ap- 
propriate if the critical loads for global buckling and local buckling both have the same order of mag- 
nitude. It should be noted that when a clamped-end simulation is used, the global buckling load is 
approximately four times the buckling load computed by VIPASA for mode (1,1). 

The modes suggested for Local Postbuckling are intended to preserve the basic shape of the 
buckling mode while allowing refinement of the shape with increasing loading. The level- 1 modes 
suggested for Local/Global Mode Interaction are intended to model the basic mechanism leading to 
imperfection sensitivity and structural collapse. The level-2 modes for Local/Global Mode Interaction 
are intended to simulate "amplitude modulation," which is the modulation (along the length of the 
structure) of the amplitude of the local-mode deflections due to the variation of bending curvature (along 
the length) due to the global-mode displacements. The strategy for selecting mode sets for local/global 
mode interaction is discussed in greater detail in [23]. 

The following additional comments apply to Table 5: 

1. It is assumed that a global mode of a symmetric structure is symmetric. 

2. If nv is even and Local/Global Mode Interaction is to be simulated, it is advised to set m, to the next 
lower (odd) integer. This is because the large bending curvatures at the mid-length of the structure 
tend to cause local-mode displacements to be maximized at the mid-length. 

3. For Local/Global Mode Interaction problems with symmetric structures, two different mode sets 
(labeled A and B) are provided with the intention that each mode set be used independently in 
separate analyses. One set models symmetric local-mode displacements, and the other models un- 
symmetric local-mode displacements. These recommendations are based on results published in 
[23], in which the direction of the global-mode response determined whether the local-mode re- 
sponse was symmetric or unsymmetric. 
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Thermal loading, as modelled with NLPAN, tends to induce in-plane normal loads. The mode sets 
suggested in Table 5 arc suitable for use in modelling thermal loading if the unit in-plane loads used in 
generating the mode shapes are similar to the in-plane reaction loads generated during thermal loading. 
Regarding pressure loading, whether or not suitable buckling mode shapes exist (for use as global shape 
functions) depends on the specific configuration considered, and this, in turn, affects the ability of 
NLPAN to model the response to pressure loads. For example, the pressure response of simply sup- 
ported rectangular plates can be accurately modelled, whereas NLPAN does not perform well in mod- 
elling the highly three-dimensional response of a pressure-loaded stiffened panel (see Section 6.3). 
NLPAN has been used to investigate the snap phenomenon in postbuckled plates [24], in which a sec- 
ondary instability (in the postbuckled regime) initiates a sudden change in the waveform. While the 
solution strategies incorporated in NLPAN are well suited to the analysis of this type of response, the 
accurate quantitative (and qualitative) prediction of secondary instabilities requires the incorporation of 
a large number of appropriately selected modes. The proper selection of these modes is a difficult 
process (see, for example, [19]) and therefore no general strategies for making such selections are offered 
here. 

53 Computer Execution Time and Memory Requirements 

The execution time for an NLPAN run is approximately proportional to the complexity of the cross 
section (in terms of the number of discretization points) and approximately doubles with each buckling 
mode added as a shape function. Typical execution times for mainframe computers and mini-computers 
are a few seconds for a single-mode analysis, and a few minutes for an analysis with ten or so modes. 

Computer memory requirements also vary with the complexity of the cross section and the number 
of buckling modes used. The NLPAN code is designed to run entirely within computer memory, so the 
size of the NLPAN analysis is limited with respect to the two characteristics mentioned. NLPAN em- 
ploys a single data vector to store, in sequence, all large data arrays, and the array dimensions are set 
based on the actual requirements needed for each specific problem. Because of this feature, NLPAN can 
use all of the computer memory which it reserves. The limit to the size of a problem which can be 
analyzed is adjusted by changing a single dimension parameter in the FORTRAN source code and re- 
compiling the code. 
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6. RESULTS 


In this section, some key features of the NLP AN code are evaluated through applications to test 
problems. First, the local-postbuckling response of a stiffened composite panel with a complex cross- 
section is investigated. Second, the nonlinear response of an imperfection-sensitive thin-blade-stiffened 
panel under uniaxial loading is explored. Third, a stiffened composite panel subjected to transverse 
pressure loading is considered. Fourth, the thermally induced buckling and postbuckling response of an 
unstiffened square composite panel is modelled. Finally, conclusions from a separately reported study 
[23] of panels and columns with constrained end-rotation are summarized. 

Some of the new features discussed in this document have been applied to test problems which are 
reported elsewhere. Results which illustrate various aspects of the nonlinear solution strategies (described 
in Appendix D) are included [16], Results from an application of NLP AN to the problem of the snap 
of a rectangular plate from one buckled waveform to another are presented in [24], 

6.1 I-Stiffened Graphite/Epoxy Panel Under Axial Compression 

NLPAN was used to model an I-stiffened graphite/epoxy panel loaded in uniaxial compression. The 
configuration is one which was tested experimentally; specifically the configuration is that of test panel 
U6 of [25]. Some features of the panel and it’s buckling response are summarized in Figure 5. The 
overall dimensions are shown in Figure 5(a). The panel featured a flat skin to which four stiffeners were 
bonded. The panel ends were mounted in potting material, and were then machined flat and parallel to 
form contact surfaces for flat-end loading. The cross-section of a representative stiffener is shown in 
Figure 5(b). A complete description of the panel is given in [25]. 

The stiffener flanges of the test panel were tapered, as depicted in Fig. 5(b). In the NLPAN analysis, 
the tapered flanges were approximated (on each side of the stiffener web) as three-step flanges by using 
three plate strips of different thicknesses. A unit-stiffener-cell representation of the panel was used, with 
symmetry conditions imposed on the skin at the edges of the unit cell. The profile of the primary 
buckling mode (as computed by VIPASA) is plotted in Figure 5(c). (The three-step flange model can 
be seen in the figure.) The buckling mode has five longitudinal halfwaves, as indicated in Figure 5(a). 

The theoretical buckling load for the full panel is reported in [25] to be 156 KN, determined using 
PASCO [13]. Reference [25] states that the mean lamina thickness was 0.014 cm The use of this mean 
thickness in the current investigation resulted in PASCO predictions of a buckling load of 210 KN. (To 
obtained this value, the critical value of axial strain computed by PASCO for the unit-cell model was 
imposed on the full-panel model assuming linear response.) The discrepancy between the two computed 
buckling loads was judged to be too large to be due to a difference in the axial buckling strain for the 
unit-cell model and the full-panel model. A second analysis was performed assuming a mean lamina 
thickness of 0.0127 cm (0.0050 in.), and this resulted in a predicted buckling load of 157 KN, almost 
exactly the value reported in [25]. Because of this agreement, the lamina thickness 0.0127 cm was used 
for the NLPAN nonlinear analysis. (Critical values of end displacement corresponding to assumed lamina 
thicknesses of 0.014 cm and 0.0127 cm were computed to be 0.094 cm and 0.077 cm, respectively, 
compared to the experimentally measured value 0.08 cm, reported to one significant digit [25], This 
result further supports the use of the smaller lamina thickness value.) 

As reported in [25], the critical load for global buckling (a single longitudinal halfwave) was well 
above the critical load for the local-buckling mode depicted in Fig 5(c), so it was assumed that the 
postbuckling response would be limited to a local -buckling type of deflection pattern. Thus, only four 
modes were incorporated as shape functions, namely (using the notation of Table 5) modes (5,1), (5,3), 
(15,1), and (15,3). The last three modes serve to refine the general shape of the first mode as the load 
increases beyond the buckling load. The unsupported length of 72.5 cm was used in the NLPAN cal- 
culations. A shape imperfection was simulated in the analysis, having the shape of the primary buckling 
mode and an amplitude of one percent of the skin thickness. All loading, displacement, and strain results 
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presented here are normalized by the theoretical values at the critical buckling load, as computed bv 
PASCO. 

Experimental and analytical results are presented in Figures 6 and 7. Results of end load versus 
end shortening are plotted in Figure 6(a), where both measures are normalized by the critical values for 
theoretical buckling. Panel failure occurred at a postbuckling load factor of 2.96 [25], The theoretical 
end-shortening values plotted are the values computed by NLPAN plus a correction for the axial 
compressibility of the potted ends, assuming that the axial strain in the ends is proportional to the axial 
load as determined by the axial stiffness before panel buckling. The NLPAN results show slightly less 
axial stiffness beyond the buckling point than was measured experimentally; this may be due to the 
difference between the clamped condition of the skin at the ends of the test panel and the simply sup- 
ported condition of the skin at the ends in the analytical model. The distribution of the longitudinal 
membrane strains in the skin across the center skin bay at the mid-length of the panel is plotted for three 
load levels in Figure 6(b). It can be seen that the NLPAN results are in good agreement with the ex- 
perimentally obtained values for all three load levels, the only appreciable disagreement being near the 
center of the bay for the higher load levels. 

Results for the variation of longitudinal surface strains with end load are plotted in Fig. 7 (with all 
values normalized). The panel locations where strains are measured are indicated in Fig. 7(a). The sur- 
face strains on the skin at the center of the panel Oocations A and B) are plotted in Figure 7(b). There 
is minor disagreement between the analysis and the experiment, but overall agreement is good. The 
opposing surface strains on and under a stiffener flange (locations C and D) are plotted in Figure 7(c). 
The experimental and analytical strain values differ by a uniform percentage over the entire load range. 
This discrepancy, which is present even in the early prebuckling regime, remains unexplained. If one 
or the other results set is scaled so that the prebuckling slopes match, then the two sets of results are in 
very close agreement. Whatever the cause of the inconsistency noted here, the agreement between the 
two results sets is still fairly good in this region of complex cross-sectional detail. 

62 Imperfection Sensitivity in a Thin-Blade-Stiffened Isotropic Panel 

NLPAN was used to model a thin-blade-stiffened isotropic panel loaded in uniaxial compression. 
The response of the configuration is sensitive to imperfections because of the interaction of the local 
and global buckling deformations. The configuration is one which was tested experimentally by 
Thompson and associates [26]. The cross-sectional proportions of a unit-stiffener-cell of the panel are 
shown in Figure 8(a). This unit-cell representation was used in the analysis; the test panel had nine skin 
bays and ten stiffeners. The panel was fabricated from epoxy resin. 

This configuration was modelled with NLPAN previously, as reported in [7]. In the previous in- 
vestigation, NLPAN was found to successfully predict imperfection sensitivity, but gave unconservative 
predictions for the limit loads of imperfect panels compared to experimental measurements. The purpose 
of revisiting the problem here is to investigate the influence of two new factors in the analysis on the 
accuracy of the predictions. The first factor is the use of the procedures described in Section 4.3.4 for 
enforcing, approximately, orthogonality between the second-order displacement fields and the buckling 
mode shapes. The second factor is the use of a mode selection strategy which enables the modelling of 
the amplitude modulation phenomenon. In this mode-selection strategy, once the local-buckling mode(s) 
to be used in the analysis is (are) identified, having a longitudinal halfwave number m, x , then additional 
local-buckling modes are incorporated which have transverse profiles similar to that (those) of the pri- 
mary local-buckling mode(s), but having longitudinal halfwave numbers (m, x - 2) and ( nu x + 2). This 
strategy is reflected in the mode-selection guidelines of Table 5, and is discussed further in [23], It was 
hoped that the presence of these two new factors would improve the agreement of the analytical results 
with the experimental data. 

Global (Euler-buckling mode) response was modelled using mode (1,1). The critical local-buckling 
mode is mode (7,1). The ratio of the critical load for local buckling P L to the critical load for Euler 
buckling P E was P JP E = 1.05. The local-buckling mode (7,3) was also deemed important so that two 
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possible modes of collapse initiation could be activated (skin buckling or stiffener buckling). Using the 
guidelines from Table 5, the following modes were selected for modelling local-buckling deformations' 

(5.1) , (5,3), (7,1), (7,3), (9,1), and (9,3). 

Local-mode imperfections of two different amplitudes were used. These were in the shape of mode 

(7.1) , with amplitudes of 2% and 10% of the skin thickness. A range of amplitudes of Euler-mode 
imperfections were used. Positive Euler-mode deflections increase the compression of the skin, and 
negative Euler-mode imperfections increase the compression of the stiffener blade. The results for limit 
load versus Euler-mode imperfection amplitude are plotted in Figure 8(b), where the limit loads are 
normalized by the theoretical critical load for Euler-mode buckling. The solid lines show the baseline 
analytical results. A second set of analytical results was generated without performing the the 
orthogonalization of the mixed-second-order displacement fields (see Section 4.3.4). These are plotted 
in Figure 8(b) with dashed lines. 

For the baseline analytical results, the limit loads for negative values of Euler-mode imperfections 
are slightly lower than those reported in [7] (not shown here), but are still unconservative. For positive 
values of Euler-mode imperfections, the baseline limit loads are actually higher (more unconservative) 
than the analytical results reported in [7], With the orthogonality condition not imposed, the predicted 
limit loads for positive Euler-mode imperfections drop sharply; they match the experimental data more 
closely for lower-amplitude Euler-mode imperfections, but diverge from the experimentally observed 
trends for larger amplitude imperfections. For negative Euler-mode imperfections, dropping the 
orthogonality condition resulted in a slight increase in the predicted limit loads. 

It was hoped that the analytical features added in the current investigation relative to the analyses 
reported in [7] would provide improved agreement with the experimental measurements. The mode- 
selection strategy used can be judged as an improvement, based simply on the argument that it enables 
the analysis to simulate the amplitude modulation which is known to occur in reality. However the im- 
position of the orthogonality condition has a mixed influence on the agreement between the analytical 
and the experimental results. Therefore, despite any theoretical justification for imposing orthogonality, 
it s not clear whether or not the accuracy of the method is improved by the practice. There do remain 
questions about some aspects of the experimental results [7]. These questions include the validity of 
assuming linear elastic material properties, and the exact shape of the local-mode imperfections, where 
the latter question concerns the fact that the imperfection amplitudes were measured only on the skin. 
Sridharan and Peng published results [4] showing good agreement between analytical results and ex- 
periment for this problem, but to obtain the results for negative Euler-mode imperfections, they used 
local-mode imperfections which highly amplify the stiffener waviness (compared to the nominal 
imperfection amplitudes) for the case of negative Euler-mode imperfections. This suggests that in the 
test panels, the stiffener waviness may have been greater than the skin waviness. In order to assess the 
accuracy of the NLPAN predictions with more certainty, it is suggested that finite-element analyses be 
performed which duplicate the NLPAN configuration and boundary conditions. This would allow an 
assessment of the NLPAN analysis without the uncertainty which accompanies experimental results. 

The use of an alternate method for generating VIPASA buckling mode shapes might improve the 
analytical predictions of NLPAN. It is noted in the discussion in [7] that a linear combination of the two 
buckling modes (7,1) and (7,3) can be used to approximately simulate isolated skin buckling or isolated 
stiffener buckling. However, the word "approximate" is important here. Mode (7,3) features large 
stiffener rolling displacements, but also includes short wavelength curvature of the skin in the transverse 
direction which would be expected to have a high level of associated strain energy. This may lead to 
suppression of the (7,3) mode, thus inhibiting the ability of the two local modes to represent two isolated 
forms of local deformation. If the two local-buckling mode shapes each represented an isolated mode 
of displacement (skin buckling or stiffener buckling) without the superfluous waviness present in the 
(7,3) mode, then the local-buckling displacements might be more easily excited, resulting in increased 
modal interaction and imperfection sensitivity. Such alternate local-buckling modes could be generated 
by exploiting the ability of PASCO to modify the transverse distribution of pre-buckling stresses so as 
to simulate transverse pressure, load eccentricity, or bowing imperfections. The local-buckling modes 
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generated in this way will tend to exhibit the types of isolated response deemed desirable here. This 
procedure has not yet been incorporated in NLPAN. 

6.3 Stiffened Composite Panel Under Pressure Loading 

NLPAN was used to model a stiffened AS4/3502 graphite/epoxy composite panel restrained at the 
edges and subjected to transverse pressure loading. The panel was rectangular with a single T-stiffener 
bonded to the center of the panel parallel to one axis. The configuration is one which was tested ex- 
perimentally, specifically Panel A of Ref. [27]. The nominal configuration of the panel test section is 
shown in Figure 9(a), and the stiffener cross-section and the laminate stacking sequences are shown in 
Figure 9(b). (The origin of the coordinate system appearing in the figure corresponds to that used in [27], 
which differs from that used in the NLPAN model.) The lamina elastic properties used in the analysis 
are also listed in Figure 9(b). The side (y-normal) edges of the panel were clamped and fixed with re- 
spect to in-plane displacements. The physical panel extended beyond the test section at each longitudinal 
end. The panel ends were not clamped at the ends of the test section; instead, the panel was supported 
against out-of-plane deflections, and the clamped condition was simulated by loading the panel with 
pressure on both sides of the end-supports. The physical ends of the panel were restrained against in- 
plane displacements. Because the length of the test section was greater than the length of the outer 
pressure-loaded bays, the effective boundary condition at the ends of the panel may have differed 
somewhat from an ideal clamped condition. 

The direction of the pressure loading is indicated in Figure 9(a). The observed panel response was 
fairly symmetrical with respect to the stiffener, so for the NLPAN analysis, only symmetric mode shapes 
were incorporated. The unit system of generalized in-plane loads used for generating the VIPASA 
buckling mode shapes consisted of a unit axial load imposed at the x-normal ends with v held to zero 
at the y-normal edges. The critical axial load for these boundary conditions was 1983 lbs. for the first 
unsymmetric mode, mode (1,1). The critical loads for symmetric modes (1,2) and (3,2) were 3265 lbs. 
and 2496 lbs., respectively. Ten modes were incorporated in the analysis, namely the first five symmetric 
modes with one longitudinal halfwave, and the first five symmetric modes with three longitudinal 
halfwaves. 

NLPAN was run first with a clamped-end condition simulated by imposing axial displacement 
constraints at the top and bottom of the stiffener blade at each end of the panel. Because the displace- 
ments computed using this representation were somewhat strange (discussed below), additional NLPAN 
runs were made, first modelling simply supported ends, then applying rotationally elastic support to the 
ends of the stiffener. The results are summarized in plots of displacement profiles presented in Figure 
10. In the figure, the rotational spring constant is denoted K„ and the normalized spring constant, defined 
in the figure, is denoted K. The transverse displacements w (positive for skin-side-out deflections) are 
normalized by the skin thickness h = .04 in.. The distribution of displacements across the width of the 
panel at the mid-length are plotted in Figure 10(a). The distribution of displacements along the length 
of the panel are plotted in Figure 10(b) for the panel centerline (under the stiffener), and in Figure 10(c) 
for the skin at y/B=0.25 . (The data plotted in Figure 6 of [27] was rescaled for plotting here in Figure 
10(b-c), to provide consistency with Figure 5 of [27], The latter figure has the correct scale; this was 
learned through an inquiry to author M.W. Hyer.) For the clamped-end simulation, the analysis predicts 
the transverse displacement of the stiffener to be essentially zero (actually, slightly negative) along the 
length of the panel, which is inconsistent with the measured response. By varying the degree of end 
constraint, the computed displacements can be brought into the ballpark of the measured displacements, 
but clearly the analysis has some shortcomings. 

The chief shortcoming identified by the author is the poor suitability of the buckling mode shapes 
of the panel for representing the pressure response. The buckling mode shapes are skin-dominated, and 
for all 5 of the modes having 3 longitudinal halfwaves, the stiffener blade remains essentially 
undisplaced compared to the skin. This causes the clamped-end boundary condition (imposed on the 
stiffener blade) to suppress the single-halfwave displacement contributions at the stiffener. These results 
expose an inherent shortcoming in the approach of NLPAN for modelling pressure-type response. The 
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deformations in the panel considered here are highly three-dimensional, and the VIPASA buckling 
modes for the panel are found to be poorly suited for representing these deformations. The pressure 
simulation has been found to work well for simple rectangular plates, and can be expected to work well 
for driving bowing-type deformation of a wide panel, because, in these cases, the buckling mode shapes 
are similar to pressure-induced displacements. 

In Figure 10(c), it can be seen that the displacement w predicted by NLP AN for the skin at the ends 
of the panel is non-zero, despite the nominal boundary condition w = 0 at the ends. These non-zero 
displacements are due to the second-order displacement fields. Despite the use of direct suppression 
of displacements as described in Sections 4. 3. 3-4. 3.4, there is a significant violation of the boundary 
condition. This occurred because the displacements at the ends are suppressed only at the node lines, 
and there is a wide expanse of skin (on each side of the stiffener) between the node line at the edge of 
the stiffener flange and the node line at the edge of the panel. A second NLP AN model was generated 
which had an additional node line in the skin on each side of the stiffener. With the displacements 
suppressed at these two additional points, the violation of the boundary condition was greatly reduced, 
and the overall displacement levels were somewhat reduced for the cases where the stiffener was not 
clamped. However there was no significant overall improvement in the analytical predictions. 

6.4 Thermally Loaded Unstiffened Composite Panel 

NLP AN was used to model the buckling and postbuckling behavior of a square, unstiffened 
graphite/epoxy panel subjected to thermal loading. The configuration is one for which analytical results 
were generated by Meyers and Hyer [28]. The plate is an eight-ply laminate with the edges simply 
supported, but with edge-normal displacements constrained to zero. The laminate stacking sequence is 
[+45/-45A)/0]s. The configuration details and material properties used in the analysis are presented in 
Figure 11(a). The plate is subjected to a uniform (change in) temperature. Meyers predicts a critical 
buckling temperature of 69.4 deg. F, whereas NLPAN predicts a buckling temperature of 7 1 .4 deg. F. 
The slight difference may be due to the fact that Meyers accounts for the laminate stiffness constants 
D 16 and whereas these values are assumed to be zero in NLPAN. 

NLPAN postbuckling analyses were performed using four different mode sets. All buckling modes 
used are sinusoidal in both the x- and y-directions. Let (m,n) denote the buckling mode with m and n 
halfwaves in the x- and y-directions, respectively. Four different mode sets were used with NLPAN, as 
listed here: 

i. 1 Mode: (1,1) 

ii. 3 Modes: (1,1), (1,3), (3,1) 

iii. 6 Modes: (1,1), (1,3), (1,5), (3,1), (3,3), (5,1) 

iv. 10 Modes: (1,1), (1,3), (1,5), (1,7), (3,1), (3.3), (3,5), (5,1), (5,3), (7,1) 

The normalized center deflection of the plate is plotted versus the normalized temperature in Figure 
11(b). All four mode sets used with NLPAN produce similar results up to a normalized temperature of 
about 1.8 . Beyond that temperature, the NLPAN results for 3, 6, and 10 modes diverge from the results 
for 1 mode, with the results for 6 and 10 modes being practically identical. The multiple-mode NLPAN 
analyses predict that the center deflection increases with temperature up to a normalized temperature of 
about 4.5, beyond which the center deflection decreases. 

The results from [28] agree with the NLPAN results up to a normalized temperature of about 1.8, 
beyond which the center deflection predicted in [28] falls progressively below that predicted by NLPAK 
The method of analysis used in [28] is more general than the method of NLPAN (for simple rectangular 
plates) in terms of modelling flexibility, but both methods share the characteristic of modelling transverse 
out-of-plane displacements using double sine functions. This author believes that the current results are 
more accurate than the results of [28] (for this specific configuration) for the following reason The as- 
sumed form for displacements used in NLPAN guarantees that for simple rectangular plates, the in-plane 
equilibrium equations are satisfied exactly for any arbitrary set of buckling modes used. This is because 
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the in-plane displacements are second-order in tenms of the modal amplitudes, so the selection of 
buckling modes determines the terms used in the in-plane displacements. With the method of [28], the 
shape functions used for in-plane displacements are selected independently of the shape functions 
(buckling modes) used for out-of-plane displacements, so that the in-plane equilibrium equations are not 
satisfied exactly unless the proper set of of terms for the in-plane displacements have been included. 
The characteristic halfwave numbers for the important in-plane displacement terms are derived from the 
sums and differences of the halfwave numbers for the buckling modes (see equations (2.7.4) and (2.7.5)), 
so that the important in-plane displacement terms do not form a contiguous group when they are ordered 
based on their characteristic halfwave numbers. Because of this, a convergence study performed by 
progressively increasing the number of in-plane displacement terms may exhibit a false convergence 
before important terms have been included. Without knowing exactly which shape functions were used 
to generate the results reported in [28], no final verdict can be reached, but the convergence of the 
NLP AN results with increasing numbers of mode shapes gives some confidence in the latter results. 

6.5 Panels and Columns with Constrained End-Rotation 

The modified-end-support modelling features described in Sections 3.3 - 3.5 were given an initial 
assessment through applications to several test problems which are described and reported in [23]. The 
conclusions noted in [23] are summarized here. 

The buckling of a slender clamped-end column was simulated. As additional displacement terms 
are incorporated into the solution procedure, the predicted buckling load converges to the theoretical 
value predicted by column theory . The NLP AN predictions of buckling load converge from above, and 
therefore, using a small number of displacement terms as is typically the practice, the results are un- 
conservative. However, for two different tests of an imperfect axially compressed T-stiffened composite 
panel with clamped ends, NLPAN predicts, with relatively good accuracy, both the mechanisms of 
structural collapse, and the limit loads suggested by the test data. A greater variety of clamped-end 
configurations need to be modelled using NLPAN in order to more fully assess the performance of the 
analytical approach. 

The mode-selection strategies discussed in [23] (and Section 5.2.2) were used in the analysis of 
axially compressed stiffened panels which were expected, based on their proportions, to exhibit local- 
global mode interaction. Amplitude modulation of the local-buckling modes during mode interaction 
was successfully modelled. Expected amplitude-modulation trends were observed for both a clamped-end 
panel and a panel with simply supported ends. 
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7. CONCLUDING REMARKS AND RECOMMENDATIONS 
7.1 Concluding Remarks 

A number of improvements and additions to the method of NLPAN were developed as a part of 
the current effort. The primary additions to the analytical capabilities are listed here: 

1. Transverse pressure loading can be modelled. 

2. Thermal loading (constant through-the thickness) can be modelled. 

3. Clamped ends, rotationally elastic end support, and eccentric end loading can be modelled. 

4. Advanced solution strategies have been implemented which allow equilibrium solution paths to be 

followed past limit points and past solution branch points of multiplicity one or two. 

5. The strategy for modelling biaxial in-plane load application has been improved, including the cor- 
rection of errors present in the original method. 

The method of NLPAN is asymptotic in nature so that solutions must be regarded as having po- 
tentially significant errors. That is not to say that accurate solutions can not be obtained in the signif- 
icantly nonlinear regime of response; the method incorporates second order contributions to 
displacements and fourth order contributions to total potential energy which are sometimes ignored by 
investigators when applying asymptotic approaches of the type used here. For simple rectangular plates, 
the NLPAN analysis degenerates to an exact series solution of the von Karman nonlinear plate equations 
(assuming that sufficiently fine discretization of they-domain is used for the numerical portions of the 
analysis). For general configurations, errors in computed solutions may be present due to the following 
factors: i) approximations made in deriving the strain-displacement relations, ii) the neglecting of dis- 
placement contributions beyond order two and energy contributions beyond order four, iii) the use of 
an inadequate number, or a poor selection, of VIPASA buckling mode shapes for use in the nonlinear 
analysis, iv) uncertainties in the method used to compute the second-order displacement fields, and v) 
the numerical error associated with the finite difference solution of the second-order displacement fields 
and the numerical integration of various functions over the transverse domain of the structure. 

For structures in which the buckling and postbuckling response is limited to local-buckling dis- 
placement shapes (little or no global-mode displacements), NLPAN gives good predictions up to loads 
of several times the buckling load. Because of its the modelling flexibility, NLPAN is well suited for 
analyzing relatively complicated cross sections for this type of response. Local/global mode interaction 
and the associated imperfection sensitivity are successfully predicted by NLPAN, although some 
questions remain about the quantitative accuracy of the predictions. Modifications were made to the 
theoretical approach, and improved mode-selection strategies were established, in attempts to improve 
the accuracy of predictions for this type of response, but the results of these effort are inconclusive. 

When NLPAN is used to model the response of a structure (panel) to transverse pressure, it has 
been found that the accuracy of the predictions is limited because of the inability of the buckling mode 
shapes to represent some types of pressure response. For simply supported rectangular plates, the 
buckling modes are well suited for modelling pressure response. The buckling modes are similarly well 
suited for modelling the pressure response of panels which buckle in a wide-column mode. However, 
for a relatively short-length, tall-stiffener panel with clamped ends and clamped edges, the buckling 
mode shapes were found to be poorly suited for modelling the highly three-dimensional response 
produced by pressure loading. 

NLPAN predictions of the post-thermal-buckling response of an unstiffened rectangular composite 
panel with constrained edges agree with other published analytical results in the early postbuckling re- 
gime. Arguments are made which support the accuracy of the NLPAN results over the other published 
results in the deeper postbuckling regime. Because the constant through-the-thickness thermal loading 
used in NLPAN acts like a form of in-plane loading (which secondarily induces bending and buckling 
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displacements), NLP AN should perform equally well in analyzing the response to either thermal loading 
or in-plane loading. 

The clamped-end modelling feature was found to give relatively accurate predictions of the response 
of a stiffened composite panel which was tested in uniaxial compression. A greater variety of test cases 
must be explored in order to more fully assess the accuracy of this modelling option. 

The advanced nonlinear solution strategies are found to generally work well, although the per- 
formance is somewhat dependent on several control parameter values and on the amplitude of 
imperfection shapes used. Sometimes numerical or approximation errors have the same effect as ge- 
ometric imperfections, and can thus cause unexpected results. Extremely small modal imperfection 
amplitudes should be avoided, because this results in solution paths with zones of extremely high cur- 
vature which the solution procedure has trouble characterizing. In general, the use of significant 
modal-imperfection amplitudes for the dominant buckling modes results in robust performance of the 
solution procedures. 

7.2 Recommendations for Future Work 

Suggestions are offered here for future work toward improving the NLP AN analysis program. Im- 
provements are sought primarily for the accuracy of the predictions for local/global mode interaction. 

1. Investigate the use of VIPASA local-buckling mode shapes which have been generated using posi- 
tive and negative eccentricities in the PASCO analysis. 

2. Compute second-order displacement fields by rigorously imposing orthogonality with respect to the 
buckling modes, and compare the fields with those computed using the methods discussed in this 
document and with results from finite element analysis. 

3. Investigate modifications to the strain-displacement relationships which may be warranted based on 
rotation and mid-surface-curvature amplitudes typically encountered. 

4. Improve the integration of the PASCO and NLPAN computer programs so as to eliminate the need 
for the redundant specification of some input parameters. 
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APPENDIX A: FORMULAE FOR THE LINEAR, UNBUCKLED SOLUTIONS 

Formulae for computing the linear, unbuckled solutions corresponding to both the unit in-plane 
loads and the unit thermal loads are discussed in this appendix. The two sections correspond to the two 
distinct load systems. 


A.l Response to the Unit In-Plane Loading 

The specific formulae for determining the solution {14} associated with the unit in-plane load sys- 
tem (see Section 2.5) are presented here. The formulae are compatible with the equations used in PASCO 
[13] but are redeveloped here for completeness in the documentation, using notation consistent with the 
present development. Biaxial loading is permitted only if there is a continuous planar skin connecting 
the boundary node lines; otherwise only uniaxial loading is permitted. Nonetheless, the solution is 
developed here assuming that a planar skin exists, and that biaxial loading is imposed, because this 
provides a solution that is applicable to all models, so long as the unit in-plane load system adheres to 
the limitations specified in Section 2.3. 

A vector {£} of length P is used to specify which plates are part of the panel skin, where P is the 
number of plate strips in the model. The elements kp of vector {it} are defined such that 

l. _ / 1 if plate strip p is part of the panel skin , , - D . , . , , . 

Kp ~ 10 otherwise (p=l,2, ... ,P) (AAA) 


The unit global load N^gl is the mean unit axial load in the longitudinal direction per unit width of 
the panel, and can be expressed as 

p 

w.1.2) 

p= 1 

where B is the reference width of the panel, N*l is the value of N x (on plate strip p) corresponding to the 
unit solution, and b is the width of the plate strip. The unit global load Nya. is the unit edge-normal load 
per unit length of the panel, and it acts on the panel skin at the boundary node lines. This load is carried 
by all plate strips in the panel skin, so that the unit y-normal stress resultant in plate strip p is given by 

(Ny L ) p = yv ^ (p = 1, 2, ... , P) (A A .3) 


The unit normal stress resultants within each plate strip are related to the unit normal mid-surface strains 
of the plate strip through the plate constitutive equations (equations (2.1.11) ): 

Wxjp = + A i2^y L )p (P= 1.2 P) (A.l .4) 

(Ny L )p = (A 12 e^ + A 22 t yi )p (p = 1 , 2 P) (A. 1 .5) 


where the unit longitudinal strain, e*., is uniform throughout the panel. 


Using equation (3) and the equation (5), the unit transverse (in-plane) strain in each plate strip can 
be expressed in terms of the unit longitudinal strain and the unit load Njcl'- 


( £ ji)p 


N yG L k P ~ £ ^(^i2)p 
(A 2l)p 


0=1,2 P) 


(A. 1.6) 


The unit longitudinal strain associated with the specified global load components can be determined by 
using equation (4) and equation (6) in equation (2) to obtain the following expression: 


^L~ 


"* L B-"yG L h 


(A. 1.7) 


where Si and s 2 are given by 
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(A 1.8) 


p 

5\ = £(*0, 

P-1 


P 


S 2 = £y«0, 

P=1 


where 


2 

C p = (^11~ A \l! A 22>p Up = ( A l2 /A 22>p (A. 1 .9) 

With the unit axial strain now known, equation (6) is used to determine the unit y-normal strain within 
each plate strip. Equations (4) and (5) are then applied to determine the unit normal in-plane stress re- 
sultants within each plate strip. 

The change in width Av t (the change in dimension between the two boundary node lines) is simply 
the sum of the changes in width of the plate strips comprising the panel skin, and this can be expressed 
mathematically as 

p 

Av L = Y, k p (U y)p (A. uo) 

P =i 

Define the mean y-normal strain of the panel skin to be 
Av l 

e yGL = ~B~ (A. l.ii) 

Using equations (6) and (10) to re-express equation (11), the mean y-normal strain in the skin can be 
expressed in terms of known values: 

*yGL = T WyG L s 3 “ V2> (A. 1 . 1 2) 

where constant is given by 

p 

*3 = '£ l k p( b/A 22)p (A. 1.1 3) 

P=1 


A few additional relationships are required for use when the boundary conditions are specified in 
particular ways. By substituting equation (7) into equation (12), the normal unit load on the side 
boundaries can be expressed as 



From equation (7) it can be determined that 


N yG L = N xGl B-Z XL ^ (A. 1.15) 

Equating the above two expressions for N^ L , the following expression for can be obtained: 

N ^l = ^- SIS3 b s ( ' 2) (A. 1.16) 

For configurations having no continuous planar skin ({*} = {0}), equation (16) degenerates to 

= (A. 1.17) 
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There are four different cases identified for sets of parameters which may be used to specify 
boundary conditions for the unbuckled panel. These cases are discussed individually below. 

Case 1) Nm. and N& specified. For this case, equations (3) through (13) provide the solution. The 
sequence of application of the equations is (8), (7), (6), (4), (3), (13), and (12). 

Case 2) specified and e,ct = 0 (Av t = 0). For this case, equation (14) is used to determine the 
effective unit load Nyct. Next, the sequence given for Case (la) above will provide the complete solution. 

Case 3) e* and e^ L specified. First, equation (16) is used to determine the effective unit load 
N*cl, then equation (14) is used to compute the effective unit load N&l. Next, the sequence given for 
Case 1) above provides the complete solution. 

Case 4) and specified. This case is required for determining the unit solution (u*} used in 
modifying the second-order displacement fields (see Section 2.7). For this case, the sequence of appli- 
cation of equations is (8), (13), (12), (6), (16), (4), and (3). 

A.2 Response to the Unit Thermal Loading 

The equations used to obtain the solution [u T ] to the unit thermal loading (see Section 3.2) are 
discussed here. Let Avv be the change in width between the two boundary node lines, for the case where 
a flat, continuous skin is present (a case where bi-axial loading is admissible). Parameter Eyjr is the 
mean y-normal strain in the skin, so that 

e yGr = Av 7 /5 (A .2.1) 

where B is the reference width. N \gt is mean load per unit panel width acting normal to the panel ends, 
and Nyor is the mean edge-normal load per unit length along the boundary node lines acting in the global 
y-direction. 

The following equations relate the various parameters applying to an individual plate strip (from 
equations (3.2.4) and (3.2.7)): 


1 

£ 

II 

(A .2. 2) 


(A.2. 3) 

1^ + A nZyr 

(A.2.4) 

N yr = ■ 4 12 e ^. + ^22^. 

(A.2.5) 

A 

where T is the unit thermal loading on the strip, other symbols are defined in Section 3.2, and the fol- 
lowing equations relate the parameters of the various plate strips to the global parameters (see equations 
(A. 1.2-3), (A.1.10), (A.2.1) ): 

1 

II 

1 

(A.2.6) 

II 

*h 

§ 

{A.2.1) 

p 

£y G T = 

p= i 

(A. 2.8) 
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where p is the plate-strip index number, P is the total number of plate strips, b is the width of an indi- 
vidual plate strip, and ^ is defined in equation (A. 1.1). 

To aid in the expression of the equations used to obtain the unit thermal solution, parameters C„, 
R p , St, S 2 , and Sj defined in Section A.l are used, along with the following additional parameters: 

p p 

^=[(^ 22 ) 0 * + 54 = '^(TbCa x ) p s 5 = ^^kp(bTF) p (A. 2.9) 

P - 1 p - 1 

The sequence of operations used to obtain the complete unit thermal response depends on which options 
are selected for control of the generalized in-plane loading (see Table 1). Four different cases are dis- 
cussed. 

Displacement control. Option 1. The homogeneous boundary conditions are given by 


6 ^ = 0 

^ = 0 (A. 2. 10) 

The following equations are applied: 

N yG T = ~S5/s 3 (A.2.11) 

= Wr l(A v) P + ( TF) p (p = 1,2 P) (A. 2. 12) 

and the remainder of the solution is obtained by application of equations (2), (3), (4), (5), and (7). 

Displacement control, Option 2. The homogeneous boundary conditions are given by 

e^=0 

N yGT = 0 (A.2.13) 

The following equations are applied: 

(ty} P = ( TF ) P (P = 1. 2, ... , P) (A. 2. 14) 

and the remainder of the solution is obtained by application of equations (2), (3), (4), (5), (7), and (8). 

Load control, Option 1 . The homogeneous boundary conditions are given by 
N*g t = 0 

ZyC T = 0 (A.2.15) 


The following equations are applied: 


„ _ S 2 S 5 + S 3 S 4 

^ 2 , 

S 2 + S X S 3 

NyG T = ( 5 4 — s \^)! s 2 

+ <” 0 , ~ R P e *r 0 = 1 ’ 2 - - • P > 

and the remainder of the solution is obtained by application of equations (2), (3), (4), and (5). 

Load control. Option 2. The homogeneous boundary conditions are given by 

^ = 0 

N yCr = 0 


(A.2.16) 

(A.2.17) 

(A.2.18) 


(A.2.19) 
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04.2.20) 


The following equations are applied: 

(ty J ) p = (f F )p- R p Ex T (p= 1,2 P ) 04.2.21) 

and the remainder of the solution is obtained by application of equations (2), (3), (4), (5), and (8). 
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APPENDIX B: ORTHOGONALITY OF THE DISPLACEMENT SHAPE FUNCTIONS 

In this appendix, the equations governing orthogonality between various shape functions sets are 
presented. First, the orthogonality condition satisfied by the buckling eigenfunctions (referred to in the 
following as buckling "modes") is developed. Second, the condition needed to enforce orthogonality 
between the buckling eigenfunctions and the second-order displacement fields is established. 


B.l Orthogonality of the Buckling Modes 


The left-hand sides of the buckling equations (equations (2.6.5), (2.6.1b), and (2.6.1c) ) are weighted 
by buckling mode components «>, v>, and vy,, respectively, and integrated over the domain of the structure. 
Because each weighted expression is uniformly zero, the integral the of weighted expressions must also 
be zero. The following equation is thus obtained: 


p . 

( J { ^xyfy) + v j(Nxyfx ^yfy t t V i«aot) + 

n= 1 N X 


' + My.jyy "f ^ I' XX ^ ^ yy^ ^ ^ ^ ^ ~ 0 


1 = 1.2, ... 
j= 1.2, ... 


(Bl) 


Equation (Bl) is manipulated by applying Green's Theorem, and invoking the definitions of { e, } , 
and {k,}, of equation (2.6.3), and {A/,}, and {A/,} of equation (2.6.2). The following equation can be 
obtained: 


P 

X V j { + + Kl N x L ( V '>x V r* + w i<x w rx> + N y^y w py^ ) 

P= 1 ^ A 

+ ^ «*[ NxUj + (Nxy. + Wx'Vi^Vj + (A 1 VX + 2 M^y + XjNxWi^Wj - M x yv j , z '] \ x = Q ^dy 

+ f + N X V j + + Myxy + XjNyW^Wj — M y Wj,y~\ = =0 


(52) 


The functional form of the buckling modes guarantees that the quantities N„, y,, w„ and M„ are all 
identically zero at x = 0 and x = L, so the second integral of equation (B2) is zero. The third integral in 
equation (B2) can be recognized to be the integral along the length of the structure of the components 
of {/.} (see equation (2.6.6) ) weighted by the components of {«*}. Using the local/global transformation 
relationships of equations (2.2. 1-2) and the definitions of {/*'"}, the third integral of equation (B2) can 
be expressed as 



[F*) T {Uj)dx 


(53) 


where N is the number of node-lines in the structure. At the non-boundary node-lines, the components 
of {Ff} are all zero (equation (2.3.1)). At the boundary node-lines, each component of {Ff} is either zero 
or its associated component of {£/f} is zero, since the buckling eigenfunctions satisfy the homogeneous 
form of whatever conditions have been specified along the boundary node-lines. Thus, expression (B3) 
is identically zero. 


Equation (B2) has now been reduced to the equation 


p . 

“J { + M) 7 ^} + Kl N xt( V i'x v j’x + w i>z w rx) + N y L w t’y w ry~} } =0 

P=1 V A 


(54) 
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Equation (B4) is reexpressed with modes i and j switched, and this equation is subtracted from equation 
(B4). Terms are eliminated from the resulting equation by recognizing that: 


{*/{£, •} = {N/{e ( } 
{MfiKj) = W/fic,} 


The following expression of the orthogonality of the buckling modes is thus obtained: 



(W/i + w i>x w rx) + N y L w by w j’y 



= 0 


1=1,2, ... 
y=i,2,... 


(B 6 ) 


A set of eigenfunctions can be obtained which satisfies the orthonormality condition 


SO. 1 ”' 


(v, 


1>X V J’X + H + A ’yWryWj'y 



(B 7) 


where 5 V is the Kroniker delta function, and a, is a constant which depends on how the eigenfunctions 
are normalized. Condition (B7) follows automatically from equation (B6) for eigensolution pairs which 
have different eigenvalues (X, * Xj). If X, = Xj but longitudinal halfwave numbers rru and m, are different, 
then equations (B7) is satisfied by virtue of the ^-dependence of the integrand. If X, = X, and m, = m, (for 
i ±j), an orthogonal set of buckling modes can be generated using the Grahm-Schmidt orthogonalization 
process. 


An alternate expression of the orthogonality condition (B7) is developed here. By applying Green's 
Theorem to equation (B7) and eliminating boundary terms which are known to be zero, the following 
equations is obtained: 


p 

X( ~f + W J< N x^bxx + UyWryyfidA \ + \ ’ n/iyctfW? \ dx = 5^ 

„ = 1 \ J A / Jo "i-"2 


(B 8) 


where is the unit global y-normal in-plane load, n, = ± l indicates the direction of the global edge- 
normal unit vector, and ru and n 2 are the index numbers for the two boundary node-lines. The second 
integral of equation (B8) is zero if the following condition is met: 


N yC L = 0 
-or- 

-or- 1^ = 0] (n = /i„ « 2 ) 


(B 9) 


The above condition is violated only if a panel has side-edges which are both i) unrestrained both with 
respect to out-of-plane deflection and out-of-plane rotation, and ii) subjected to y-normal in-plane load- 
ing; commonly encountered configurations generally do satisfy equation (B9). If it is assumed that the 
conditions of equation (B9) are satisfied, the orthogonality condition of equation (B7) has the following 
equivalent form: 



Xj ^'rxx) 'f ^j^x^bxx 


+ 


)]dA^ = -S i/ii 


(BIO) 


B.2 Orthogonality Between the Buckling Modes and the Second-Order Fields 

Assume that it is desired to enforce orthogonality between each second-order displacement field 
{«.;} an( i each buckling mode {«*). In the orthogonality condition for the buckling modes, equation (B7), 
replace {«,} with {«,;): 
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p / 

y, ^Xf( V ij<x V k<x + w ij>x w bx) + Ny • / } V ij<y w by~\dA^ = 0 


O’, 7 = 1,2, ...) 
(*=1,2, ...) 


By applying Green's Theorem to equation (B 11), cancelling boundary terns that are identically zero, and 
assuming that equation (B9) is satisfied, equation (B 11) can be converted into two alternate forms: 


V, ( J I- v *(^c L v i>*xx) + ^k^x^ipxx + ^yL^ij'yy) 3 ^ / = 1 

p=\ Xa / p 

p / 

~ , f J t v ij^x { ^bxx) + ^if^x^k'xx + ^k'yy)\dA ^ 

P p A /p 

+ X n * N *S v ij v *x + w tJ w bx) \ x = Ql dy^ = 0 


(512) 


(513) 


Functions v v and w v do not necessarily go to zero x = 0,L; however, according to the nominal boundary 
conditions, the transverse displacements v and w should be approximately zero at x = 0, L. Assume that 
the solution procedure used to compute { «,>} is successful in assuring that this boundary condition is 
approximately satisfied, so that the boundary term in equation (B13) can be neglected. Equation (B13) 
then becomes 


p 

^ ( j C V i/.^x L v k<xz) + W i/~^z : W k’xx + NyW k ,yy)~\dA j =0 
P=1 V/4 Sp 


Thus, equations (B 11), (B12), and (B14) are equivalent expressions for orthogonality of the functions 
[Uij] with respect to functions {«,}. 
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APPENDIX C: AN EXPANSION OF dn 3 


In terms of the notation used in equations (2.8.29-30), 5ttj can be written as 
8 % = 5qfljq k (C ijk + Xcfji) (Cl) 


for the case of a perfect structure. In order to keep the following derivation manageable, attention is 
limited to the case where displacement control is used for the in-plane loading (CONTRL='D' in Table 
1), and where coefficients B, (1=1,2, ... ) appearing in equation (2.6.12) are zero. With these re- 
strictions, equation (Cl) can be expressed in terms of "primitive" coefficients as 

8k 3 = &7,<7/7*[(2c£ + C;*) + X(cjj k + 2 c£*)] (C2) 


The primitive coefficients have the following definitions: 




/ P 



[N L ) T {t ijk }dA 



(C3) 


where [N l } is defined in equation (2.5.3), {A() and {M,} are defined in equations (2.6.2), {e, ; } and 
(K.yJ are defined in equations (2.7.3), and {e*} is defined in equation (2.8.6). 


Through applications of Green's theorem and cancellation of boundary terms at the jc-normal ends 
which are known to be zero (because of the function forms of { 14 } and {«,>)), the following two equiv- 
alent expressions are developed for Cj: 


C-j= - l<uij, L(N*)> + <u t , N(N k , «,)>] 


m 


+ t 1 

= - <u k , L(N i} )> 

P 


n X-^xy i ^ij + 

k u ij + N yk Vy (2Mxy k 'x My k ty)Wjj AfyW^y^\ 
— [NyUj,^ + (NfyWj^ + NyWj, y )W~\ | I J,= 0,6 d*) 


(C 4 ) 


M x .yv k , x 




; u k + N yVk + (ZMxy^x + M y .., y )w k - MyW k ,y~\ I v _ n h dx 


y,jk'yJ 1 j, = 0.6 


■). 


where the inner products and operators which appear are defined in equations (4.2. 1-3). Similarly, the 
following two equivalent expressions are developed for C%: 


Ci jk = - <Uj k , N(N l , «,)> 

p 


^ ( l + W, ' xWjk> ^ - <U, ^ + ^ n > N > l W ‘'> W jk\y = 0.b 


(C 5) 


= -<u if N(N l 


. «/*)> + m riyNy 


y L W jlcy W i 1 ^ = 0,b 


*). 
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Equation (C2) is evaluated using equations (C4) and (C5). For Cj, the first form in equation (C4) 
housed, and for C the second form is used. For C&, the second form in equation (C5) is used, and for 
Qk, the first form is used. Terms on the .y-normal boundaries of plate strips are transformed to refer to 
the node-lines and global coordinate directions, and the contributions to the generalized node-line force 
resultants are summed, where appropriate. The following equation is obtained: 


5*3 = 5?,<7/7*{ - <«,, L(N jk ) + N(N k , u ; ) + lN(N L , u jk )> - 2<u I> , L(/V*) + XN(N l , «*)> 

*i« xm 0 l dy^ + Xjf l ^ )r(/ >} ctc 

la t V i/^xy k + ^^x L V k’x) + w ij fMx^x + ^^xy k ’y + 3 I x « OJ, ^ ^ 

+ { + 2 (^ “ 1 -a. -z ^ } 


(C 6) 


Several simplifications are made to the above equation. First, the integrals containing {(£}, {£/*}, 
{/>}, and (FT) are deleted, because the homogeneous node-line boundary conditions specify that one 
or the other term of each product is zero. Second, while the functional form for {«,>} does not guarantee 
it to be so, it is assumed that at x = 0, L the quantities A/*,, v,>, and w,>- are approximately zero (consistent 
with the buckling mode characteristics ,■= v, = w, = 0 at x = Q, L). This eliminates the integrals at the 
x-normal ends of the structure. Third, the structural configuration is assumed to satisfy Appendix B 
equation (B9), thus eliminating the integral at the boundary node lines «, and n 2 . (This assumption states 
that the configuration does not have a free y-normal edge with y-normal in-plane loading.) Finally, the 
buckling equations as written in equation (4.2.5) are used to further simplify the second inner product 
in equation (C6). The following equation is obtained: 

5*3 = 5 MMflk { - <«,. lUN jk ) + N(N k , uj) + i- N(N jt Uk ) + XN(N l , «;*)]> 
-2(X-X ik )<u (/ ,/V(iV L ,« jk )>} 

where the expression [<7/7*N(M, My)] of equation (C6) was manipulated to obtain symmetry in indices j 
and k. 
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APPENDIX D. DESCRIPTION OF THE NONLINEAR SOLUTION STRATEGIES 

D.l. Form of Equations Governing Equilibrium 

It is assumed that the total potential energy of a structure can be expressed in the following algebraic 
form: 


JC = 7C(<7, a, p) 

-K 0 + Qi (Ai + ttA/ + P^f) + qfljiAy + aA “ + pA ( y) (0 i ) 

+ WM A ‘jk + ^ijk + P A? jk ) + qflflkq[(A ljki + a4 + p A$ a ) 

where summation over i,j, k, and / is implied, q is a vector of generalized coordinates, a and P are 
generalized load parameters, the sub- and super-scripted coefficients A are constant, and the term n„ has 
no dependence on the generalized coordinates. It is assumed that the generalized coordinates and load 
parameters in equation (Dl) have been normalized so that they take on values of order of magnitude 
unity in the course of an analysis. The method described here is not limited to total potential energy 
expressions with two load parameters and fourth order terms; these specific characteristics are adopted 
for demonstration purposes. 


The equations governing equilibrium are obtained by imposing a stationary total potential energy 
condition, expressed as 

Si(q . a, P) = 0 (i=l, 2, ... ) (02) 

where 


" dq, 

= (fl, + a B “ + psf) + q/B 0 + a£,“ + pfif) 

+ qM B ijk + a tfk + pfljt) + qflMB l]kl + aB? Jkl + p Bf jU ) 

The newly introduced coefficients appearing in the above equation are given by, for example, 

B i = A i B ijk = A ijk + A Jlk + A jki 

B ij = A ij + Aji B ijkJ = A iJkl + A ]iU + A ]kil + Aj Ui 


( D3 ) 


(D4) 


The two load parameters a and P can be controlled asynchronously using a single load parameter 
X. A series of K load ranges is specified in terms of target values for a and P: (0, 0), (a,, p,), (a 2 , p 2 ), 
... , (a*, p*). Over the k* load range, a and P vary linearly with X as X increases from 0 to 1: 


(p) ■ (tO + 


0<X< 1 


(D5) 


For the load interval, the equilibrium equations then take the form 

fi(q< X) — (C, + XC , ) + qfC l} + XCjj) + qjq k (C tjk + XC t]k ) + q/Ji^iiC^u + XC ljk j) 

= 0 0 = 1.2 M) 


(D6) 


where for example, 

Q = Bi + & k -\ B* + p A _ , Btf 

cN(a*-a*-i)fi“ + (P*-P*_,)fif 


and where it is assumed henceforth that a finite basis of M generalized coordinates is in use. 
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D.2. Arc-Length Control Method 

The implementation of the Riks-Wempner arc-length control method described here was guided in 
large part by the presentation given in [Dl]. Concepts presented in [Dl] are also used here in a proce- 
dure for locating and classifying critical stability points. 

_ D.2.1 . Arc-length parameter "s". Let ( q , X) be a known solution to equation (D6), and let 
(q + £, X + 8) be a new solution which is sought in the vicinity of the known solution. The new solution 
satisfies the equation 


ffi + Z,,X + 8) = 0 (z = l,2 M) 


(0 8 ) 


The independent parameter s is introduced, so that q = q(s) and X = X(s). Parameter s is the arc- 
length measure in ( M + l)-dimensional space for an equilibrium solution path, and is governed by the 
equation 


qjqj+X =1 (D9) 

where summation over j is implied, and where q t = dqjds and X = dX/ds. A Taylor sejies expansion about 
the starting solution is used to obtain an expression for the incremental solution (£, 8) in terms of the 
an arc -length increment A s: 


£, = qAs + OiAs 2 ) , 8 = X As + O (As 2 ) 


(DIO) 


D.2.2. Determination of the derivatives (q, X). By expressing each equation of equation (D8) as a 
Taylor series expansion about the solution (q, X), applying the substitution of equations (DIO), and taking 
the limit as As approaches zero, it can be shown that [Dl] 


D ij q j + D?X = 0 , (i = l,2 M) 

In equation (Dll), D, y is the tangent stiffness matrix, given by 


D ii = 


dqj 


Cij + Qk^Ciji + Cjy) + q^fCiju + Cm + Cjjj) 


(Dll) 


(£> 12 ) 


(?. *■) 

where, for example, 


Cijk ~ Cijk i/Jt 


and vector Df is defined by 

dX , , 

(?•>.) 


d. 8 = 


= c} + qf) + qjq k cf jk + qfl^ju 


(013) 


(014) 


The derivatives {q, X) are determined by arbitrarily specifying one of the M + 1 unknown values, solving 
for the remaining values using equation (Dll), then scaling the solution so that equation (D9) is satisfied. 
The values (q, X) thus obtained are ambiguous to the extent of a factor ±1. For a starting solution 
(A. = 0), X is positive. In a solution stepping procedure, it is assumed that the direction of the vector 
(q,X) does not change radically between successive solutions, so that at the n* equilibrium solution 
(< 7 * , X"), the sense of the vector (q*, X’) is determined by the following requirement: 


q- l qj+X n ~'X n >0 


(015) 

The use of this simple technique negates the criticism of the arc-length control method, issued in [D2], 
that the sign of X is ambiguous. 
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D.2.3. Solution stepping using arc-length control. Using equations (D9.10) it can be determined that 
for a sufficiently small increment As the following approximate relationship holds: 

As « qfa + X& (D16) 

where summation over j is implied. Arc-length control is imposed by specifying the increment As from 
the n * solution to the (n + 1)** solution, and assuming that the approximate relationship of equation (D16) 
is an exact one. The solution is then governed by both equation (D8) and equation (D16). Newton- 
Raphson iteration is used to determine the new solution in the following way. Where (q*, X r ) denotes 


the r * estimate of the new solution, the iterative solution procedure can be written as 

(q r + l , X r + ! ) = fq, X) + (Aq + \ AX + *) (DM) 

where the correction (Aq'^', AX ,+ 1 ) is obtained by solving the linear system of equations 

D ij + (D?Y AX + 1 = - R- (/=1, 2, ... , M) (D 1 8a) 

q,Aq' +l + XAX +i = As - \_qfq] - qj) + X(X - V)] (D 1 8 b) 

where matrix Dj and vector (Dfy are evaluated at (q r , X), and Ri is the residual error vector for the r* 
estimate: 

R r , =/<^. X) (D19) 

D.2.4. Stability criterion; properties and classification of critical points. The following eigenvalue 
problem is evaluated for use in assessing the stability of an equilibrium solution: 

([D] - co*[/]){<t>*) = {0} (Ic = 1, 2 M) (D20) 


where [£)] is the tangent stiffness matrix and [/] is the identity matrix. Eigenvalues to* are assumed 
to be ordered according to increasing value. When all eigenvalues are positive, the tangent stiffness 
matrix is positive definite, the total potential energy is a local minimum, and the equilibrium state is 
stable. At a critical stability point (a limit point or bifurcation point) the first eigenvalue co, is zero. 
Consider co t to be a function of the path parameter s, and let s* be the value of s at a critical point. The 


critical stability point then has the property 

«>,(■**) = 0 (D21) 

A limit point has the property 

Us*) = 0 (D22) 

whereas a bifurcation point has the property [D 1 ] : 

B,(s*) = 0 (D23) 

where B k (s) is defined as 

B k (s) = (4>*D/) L (lc= 1,2 M) (D24) 


It is noted that at a limit point, the path tangent q is a scalar multiple of the eigenvector <j>\ This can 
be seen by considering equations (D 11,20-22). 

D.2.5. Detection and classification of forward critical points. The criteria given in equations 
(D21-23) for identifying and classifying a critical point apply only at the critical point itself. A method 
is described here for detecting and classifying a "forward" critical point, meaning one which is being 
approached in the process of solution stepping. The values of o>i, X and B x are computed for each dis- 
crete equilibrium solution. In the vicinity of the most recently obtained solution, the three functions 
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t»i(s). Ms), and 61 ( 5 ) are approximated as quadratic functions of s by fitting curves to the parameter 
values from the three most recently obtained solutions. The equations C 0 i(j*) = 0, X(s,*) - 0 , and 
Bfs b *) = 0 are then solved to obtain extrapolated values for s*, s,* t and s**, which are the values of s 
at the next forward critical point, limit point, and bifurcation point, respectively. Each of the parameters 
s*. Si*, and s t * is set equal to the smallest positive real solution of its associated quadratic equation, if 
such a solution exists; otherwise the parameter is set to a large, positive real number. 

Define As*, As,*, and A s b *, to be the incremental values from the current value s to the values s*. 
Si*, and s b *, respectively. If, As,* approaches As* as a critical point is approached (ie. as As* approaches 
zero) then the critical point is classified as a limit point, and if As** approaches As*, the critical point 
is classified as a bifurcation point. (A critical point may be both a limit point and a bifurcation point) 
In order to evaluate these trends, it is necessary to compute several equilibrium solutions which approach 
the critical point without overshooting it. To do this, Riks [Dl] suggests selecting As to be As = cAs* 
(where c is a factor less than one and greater than zero) until the extrapolated value As* is smaller than 
a specified cutoff value. Because the critical point can be approached only to within some finite incre- 
ment As*, there is some uncertainty involved in evaluating the trends of the extrapolated values. The 
approach used here is to inspect the ratios (As ; */As*) and (As**/As*) evaluated at the smallest value of 
As* used. If a ratio has an order of magnitude of unity (a value less than, say, 3) then it is assumed that 
the associated values are converging. The use of quadratic interpolation functions has been found by the 
author to provide a considerably more robust method of locating critical points than when the linear 
interpolation functions suggested in [Dl] are used. 

Plate problems have been encountered in which the parameter B, is uniformly zero along a 
postbuckling equilibrium path. In this case the extrapolation procedure described for determining Ay** 
is not appropriate, and any critical point encountered is a bifurcation point. 

While equation (D21) defines a critical stability state (<Di = 0), other singular points may be en- 
countered where the first eigenvalue g>i is negative, but another eigenvalue to* is zero for some value k 
greater than one. The criteria of equations (D22) and (D23) can still be used to classify the singular point 
as a limit point and/or bifurcation point, except that B k is used in equation (D23) rather than B t . 

D.2.6. Computing the solution at a critical point. If it is determined that a critical point is not a 
bifurcation point, then the critical equilibrium solution (a limit point) is determined by selecting As of 
equation (D18b) equal to the extrapolated value of As*, once the latter measure is sufficiently small, and 
using the iterative solution procedure already described [Dl]. If it is determined that a critical point is 
a bifurcation point, arc-length control is not well suited for centering on the critical point because of the 
solution branching. It is the author’s experience that it is possible to approach a critical point quite closely 
using arc -length control, to the point that the value of As* is very small compared to unity. It thus seems 
sufficiently accurate to perform a simple extrapolation to determine the bifurcation point by using a 
first-order approximation of equations (DIO) with Ay = Ay*. 

D.3. Control of Solution Branching with Thurston's Method 

This section describes an implementation of Thurston's method [D3], used for analyzing solution 
branching in the vicinity of a bifurcation point which has been identified and located using the proce- 
dures described in Section D.2. In [D3], Thurston's method is applied to the differential equations 
governing equilibrium before any necessary discretization of the structural response has been performed. 
Here, the discretization has already been performed in obtaining equations (D 6 ) which govern equilib- 
rium, and Thurston's method is applied directly to these equations. 

D.3.1. Transformation of the equilibrium equations. Let (q, X) be a known exact or approximate 
solution to equations (D 6 ), and let (q + X + 5) be a new solution which is sought in the vicinity of the 
known solution. The latter solution is governed by equation (D 8 ). By expressing equation (D 8 ) in the 
expanded form of equation (D 6 ) and grouping terms based on their power in and 8 , the following 
equation is obtained: 
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( 025 ) 


<D, + 5 of) + %fPjj + 50?) + !&<D* + 6 Dfo + + SC?„) = 0 

(/ = 1, 2, ... , Af) 

where £>,* and D v are given in equations (D14) and (12), respectively and 

a =/<?. >•) 

D ?j = 4 + <^4* + Cj&> + wXCjn + 4, + (*/) 

<=cJ t +< ?/ (cJ u +cJ*+cJ*) 

— Cijk ■*" QliCijkl ^ijlk + ^j/yt) + ^E>ijk 

Dijkl = Q/*/ + 


(D26) 


The following eigenvalue problem is derived from the linear terms of equations (D25): 

(M + 8*[D 6 ])(6*) = {0} , k = 1, 2 M (027) 

where [o] and [D s ] are the matrices of coefficients £> v and Dg, respectively. Equation (D27) differs 
form the eigenvalue problem of equation (D20) by the presence of the matrix [ -£>*] in place of [/]. 
The eigenvalues 5* are numbered in order of ascending value. When all eigenvalues are distinct, the 
eigenvectors form an orthogonal set, and the eigenvectors are scaled to meet the following 
orthonormality condition: 


0^0* = -8^ 


(028) 


where 8,y is the Kroniker delta function. For multiple eigenvalues which coincide, the associated 
eigenvectors must be orthogonalized before normalization, using, for example, the Grahm-Schmidt 
orthogonalization procedure [D4], 


The incremental solution E, is expressed as a series in the eigenvectors, 0*: 
l = ^ k °k (029) 

where summation over k is implied, and the coefficients a* (k = 1, 2, ... , Af) are initially unknown. The 
following relationship is established through the use of equations (D27-29): 

0”(Oy + 5 D^j = - (5 - hja m (m = 1 , 2, ... , M) (030) 


where summation of repeated indices is implied except for index m. Equation (D25) is now restated, 
eliminating parameters using the substitution of equation (D29). For each value of 
m (m= 1, 2, ... , Af), the i* equation is weighted by 0T and the resulting expressions are summed over 
i. Equation (D30) is used to simplify the equations thus obtained, providing the following transformed 
equations: 

^(8 - 8 J = (E m + 84) + (£w + 8E*„) + a s (E mnrs + SE^J j ) 

(m = 1, 2, .... A/) 


where summation over n, r, and s is implied, and the new coefficients are given by 


(£ m ,4) = 0 r(O ( ,of) 

tfw eLs) = 'fi&kPyu. c} jU ) 


(E mnr ,El r ) = QTQ]Q r k (D ljh Df jk ) 


(032) 


where summation over i,j, k, and / is implied. 


D.3.2. Analysis of solution branching. The form of equations (D31) is a generalization of the form 
of equations obtained by Thurston in [D5] for the analysis of compressively loaded plates, so comments 
made in [D5] will be used to guide the use of equations (D31). Near a limit point or an isolated 
bifurcation point, one of the eigenvalues, say 8„ will be close to zero and much less in amplitude than 
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the remaining eigenvalues. In the vicinity of the known solution, both 6 and the quantity (8 - 8,) are 
small. With this in mind, and in view of the form of equations (D31), Thurston hypothesizes that a, will 
be much larger in amplitude than the remaining M — 1 coefficients a m . An approximate solution is thus 
governed by equation number t of equations (D31) with only the variable parameters a, and 8 retained. 
By specifying a value for a,, the corresponding value of the load increment 8 can be computed directly. 
Solutions for both positive and negative values of a, are generated in order to identify both solution 
branch directions. 

When two eigenvalues 8, and 8v are both very close to zero (such as near a point of approximately 
simultaneous buckling in two different mode shapes) the reasoning used above is extended to suggest 
that an approximate solution is governed by two equations, numbers t and v of equations (D31), with 
only the variable parameters a,, a „ and 8 retained. The solution procedure used here is to specify a small 
numerical value for either a, or a,, solve the two equations for 8, and then equate the two expressions 
for 8. A fifth degree polynomial equation in the remaining unknown parameter (a, or a,) is obtained. A 
polynomial root solver is used to determine all real solutions to the polynomial equation, then the cor- 
responding values of 8 are computed. Solutions are generated for both positive and negative values of 
the specified parameter (a, or a») in order to identify all solution branch directions. The analysis of si- 
multaneous bifurcation in three or more modes has not been attempted as a part of this work. 

The set of approximate solutions obtained using a one- or two-mode branching analysis are assessed 
for physical significance in order to guide the selection of a particular solution branch to follow. The 
approximate solution on the selected branch is then refined using the following procedure. The approx- 
imate post-bifurcation solution is known in terms of a value 8 (henceforth denoted 8*) and one or two 
non-zero parameters a,. Using equation (D29), the corresponding approximate solution is generated. 
The arc-length increment corresponding to (£*, 8*) is given by 

A/ = V^? + (SV (D33) 

and the approximate path derivatives along the new branch are given by 

t=?/As b , X=8 b /As b (£>34) 

The arc-length control method is used to refine the solution by taking the extrapolated critical-point 
solution (q*, X*) as the starting solution, (q, X) of equation (D34) as the path tangent, Ay* as the arc- 
length increment, and ( q * + £*, X* + 8*) as the initial guess for the new solution. 
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USER INSTRUCTIONS FOR THE NLPAN COMPUTER PROGRAM 


Introduction 

This document provides user instructions for the NLPAN code as it is currently configured under the 
Unix operating system. Detailed line-by-line instructions for preparing input are preceded by discussions 
of general aspects of running NLPAN, general aspects of specifying geometry, various operating modes 
of NLPAN, convergence considerations, and possible problems which may be encountered. Finally, some 
example input files are listed. The user must be familiar with PASCO input procedures, because a PASCO 
input file must be formed for the structure to be analyzed. PASCO input requirements are described in 
[El], and are not included in this document. 

General Aspects of Running NLPAN 

All FORTRAN source code for NLPAN resides in a single sub-directory, along with a ’makefile’ used 
to compile and link the program. These files can be copied to a user’s own account if it is desired or 
necessary to alter the program in order to increase array sizes, generate specialized output, etc.. 

The source code for PASCO has been incorporated into the NLPAN code, and a conventional PASCO 
input file is used to specify much of the configuration geometry for an NLPAN analysis. However the 
PASCO source code has been altered and reduced in its capabilities based on the needs of NLPAN, and 
thus PASCO can not be used as a stand-alone program from within NLPAN. 

NLPAN runs entirely within computer memory, with the exception that a couple of data storage 
functions use disk files. This limits the numerical size of problems which can be considered. All 
dimensions for data arrays are set in a PARAMETER statement in the file ’param.f’. A single mass 
storage array is apportioned as needed to store many large arrays without wasted space. The dimension 
of the mass storage array is set by parameter MMASS. When dimensions have been exceeded, error 
messages are printed. The dimensions can be adjusted in the PARAMETER statement of file ’param.f’, 
and then the program is recompiled using the ’make’ command with the provided ’makefile’. The standard 
output file ’nlpan.out’ includes a listing of dimension NMASS, which is the length of the mass storage 
vector required for a particular problem. 

Depending on the problem size, NLPAN execution typically takes from a few seconds to a few 
minutes to run. The execution time is approximately proportional to the complexity of the cross section 
modelled (in terms of the number of plate strips and the number of discretization intervals on each plate 
strip) and approximately doubles with each additional VIPASA buckling mode included as a shape 
function. 

NLPAN is executed interactively from within a data sub-directory containing the input files for a 
particular configuration to be modeled. The run command is ’(path)nlpan’ where ’(path)’ is the path 
specifying the location of the executable file. The Unix command ’alias’ can be used to define a simple 
input string, such as ’nlpan’, which initiates program execution without the need to type the complete path 
each time the program is run. The data sub-directory must contain three input files: 

’pasco.in’ PASCO input file 

’nlpan.inl’ Basic NLPAN input including modelling options and geometric data. 

’nlpan.in2’ Imperfection amplitudes and control parameters for the final nonlinear analysis. 

All working files and output files are written to the data sub-directory using generic names, so it is natural 
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to use a separate sub-directory for each configuration analyzed. The standard output files are: 
’pasco.out’ Output from PASCO 

’nlpan.out’ Standard output from NLP AN, including basic load-step output 
’nlpan.outr’ Output of analysis-control parameters used in solution stepping 

Output files are overwritten with each successive run, so files must be renamed to be saved. Additional 
output files are generated as selected by the user using parameters described in the section "Line-by-Line 
Input Instructions." Several files are created for passing data between different program units, including 
’stiffin’, modes.in’, ’reverse.in’. If intermediate calculations are saved for restarting, the files 
’nlpan.binstor’ and ’nlpan.binind’ are created. 

General Aspects of Specifying Geometry 

The first step is to form the PASCO input file, which is given the name ’pasco.in’. The following 
limitations are placed on the PASCO model: 

1) The configuration must be defined using both HCARD and ICARD input, with ICREP=1 and 
NOBAY=l. This suppresses substructuring in the VIPASA buckling analysis so that the buckling 
modes returned from VIPASA include information at all of the node lines. 

2) Only in-plane normal loads NX(1) and NY(1) may be specified. The PASCO features of shear 
loading, pressure loading, load eccentricity, bowing imperfection, and vibration frequency should not 
be used. 

3) PASCO input data related to design optimization procedures are not used. 

4) The following PASCO input parameters should NOT be present: MINLAM, NEIG, and NLAM. 

5) Set CONV1=50000„ FREQ=0., MAXJJJ=0, LINK=1 

In defining the PASCO geometric model, a cross sectional configuration is defined in terms of a set 
of numbered plate strips and a set of numbered node lines. An example of a model with 7 plate strips and 
6 node lines is shown in Figure El(a). Each plate strip has a local Y coordinate axis (indicated by arrows 
in Figure El (a)) originating at one of the edges of the strip. Each plate strip may be rotated by an angle 
MU relative to the global coordinate directions, and a plate edge may be offset from its associate node 
line, as specified by the values EY and EZ defined in Figure El(b). The connectivity of the model and 
the rotation angles MU and offset measures EY and EZ are not passed automatically from PASCO to 
NLPAN, and thus the information must be repeated (in different form) in input file ’nlpan.ini’. The 
connectivity is specified by using the NOD array in ’nlpaainl’, which lists, for each plate strip, the i.d. 
numbers of the node lines at Y=0 and Y=B, respectively, where B is the width of the plate strip. The array 
NOD, the angles MU, and eccentricity measures EZ for the configuration of Figure El(a) are listed in the 
figure (EY is zero for all plate strips). The details of specifying these parameters are described in the 
section "Line-by-Line Input Instructions." The NLPAN input files for the configuration of Figure El(a) 
are included in the section "Example Input Files." The following data are passed automatically from 
PASCO to NLPAN: plate-strip widths, global unit loads, laminate configurations, reference length and 
width, and material properties including thermal expansion coefficients. 
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(a) An example PASCO/NLPAN cross-sectional representation showing 
plate strip numbers (boxed), node-line numbers (circled), global and local 
coordinate directions, and geometry definition parameters. 



(b) Orientation of the edge of a plate strip relative to a node line. 


Figure El. Specification of cross-sectional geometry. 
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One important modelling detail requires special discussion. In both PASCO and NLP AN, a set of 
plate-strip properties is stored only for unique combinations of the following three parameters: 

i) wall Qaminate) properties 

ii) plate-strip width 

iii) unit in-plane loads 

In NLPAN, the user must use the ILBWAL(IP) input to assign to each plate strip the index number of 
this property set, and thus the user must anticipate the indexing used by PASCO. This is not difficult. 
Starting with plate strip # 1 and proceeding through the entire set of plate strips, assign index numbers 
ILBWAL starting with 1. If the trio of properties of a plate strip match the values encountered on a 
previous plate strip, use the previously assigned value of ILBWAL, otherwise assign the next new value 
of ILBWAL. The vector ILBWAL(IP) for the configuration of Figure El(a) is (1,2,2, 1,3, 3, 4), where the 
laminate of the stiffener flanges is different from the laminate of the skin. On occasion, PASCO does not 
recognize matching plate strips. A FORTRAN error results if there is a discrepancy between the PASCO 
indexing and the ILBWAL input. The file ’modes.in’ should be inspected. At the top is variable NLB WAL 
which is the maximum value of ILBWAL used by PASCO. Lines near the top list parameters ILBW, IW, 
B, NXL, NYL. IW is the wall (laminate) index number, ILBW is the index referred to by ILBWAL(IP), 
B is the strip width, and NXL and NYL are the unit loads. These values reflect the indexing used by 
PASCO. 

In the PASCO input file, the four boundary condition components at the node-lines are specified in 
the order dW/dY, W, V, and U. The boundary conditions must be specified again in input file ’nlpan.ini ’ 
using the BCVEC parameter. The components of BCVEC are reversed compared to PASCO: (U : Fx), 
(V : Fy), (W : Fz), and (dW/dY : M). 

NLPAN has features for modelling boundary conditions other than simple support at the longitudinal 
ends. Axial displacement constraints can be imposed using NCU/IPCU/YCU input, and slope constraints 
dW/dX can be imposed using NCW/IPCW/Y CW input. The generalized displacement constraints are 
satisfied exactly at discreet points; care must be taken not to over-constrain a problem. Rotationally elastic 
support can be modelled using either NSPU/IPSPU/YSPU/KSPU or NSPW/IPSPW/YSPW/KSPW input. 
The options provided for modelling end-support are not intended to enable exact modelling of boundary 
conditions at the longitudinal ends, but are rather intended to simulate the effects of the end support 
conditions on the global behavior of the structure. Axial displacement constraints are used in the example 
input files listed at the end of this document. 

NLPAN Operating Modes 

NLPAN has several special operating modes, and to understand them, it helps to understand the 
general stages of the NLPAN analysis: 

1. Computation of the linear, unbuckled response to the unit in-plane loads. 

2. Computation of VIPASA buckling eigensolutions, and selection of the buckling mode set for use 

in the nonlinear analysis. 

3. Computation of the shape functions for the second-order displacement fields. 

4. Calculation of the coefficients for the nonlinear algebraic equations governing equilibrium. 

5. Solution of the nonlinear algebraic equations governing equilibrium. 
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Linear solution pre-processor. NLP AN provides the ability to specify both global in-plane load 
components in terms of either load measures or displacement measures. If it is desired to use the two 
edge-normal load components (named NX(1) and NY(1) in ’pasco.in’) then PASCO will compute the 
linear pre-buckling solution, which is passed on to NLP AN. If it is desired to specify the unit end- 
shortening and/or the unit width-shortening, NLP AN is used as a pre-processor to determine the equivalent 
loads NX(1) and NY(1) required for PASCO input. This is done by setting ILNPRT=2 (pre-processor 
mode) in file ’nlpan.ini ’ and providing the appropriate values for ILSET, NXGL, NYGL, EPXL, and 
EPYLG in ’nlpaainl’. NLP AN then computes the corresponding values NX(1) and NY(1), prints these 
values to file ’nlpan.out’, and halts program execution. The provided values can then be put into ’pasco.in’ 
(using a text editor) for use by PASCO, in order to get the proper pre-buckling state for the buckling 
calculations. ILNPRT should then be set to zero for the subsequent NLPAN run. 

For multiply-connected cross-sections, the PASCO linear solution is computed without strict 
enforcement of displacement compatibility between multiple plate strips. By setting COMPAT=’Y’ in 
’nlpan.in’, the linear solution can be determined which satisfies displacement compatibility for 
configurations with co-planar plate strips, such as skins with attached flanges. First set ILNPRT=2 (pre- 
processor mode) in ’nlpan.in’. This will cause the PASCO input parameters NX(1), NY(1), and FNY(I) 
to be printed to file ’nlpan.out’, and NLPAN program execution will terminate. (If NY(1) is zero in 
’pasco.in’, then NLPAN will also specify a small non-zero value NY(1) for use in PASCO.) The values 
NX(1), NY(1), and FNY(I) are then added to ’pasco.in’ using a text editor, and ILNPRT should be set 
to zero in ’nlpan.in’ for subsequent runs. The linear solution obtained using this method may not be a true 
equilibrium solution because of the presence of unbalanced moments along the node lines. 

Review and selection of VIPASA buckling solutions. The set of VIPASA buckling modes used by 
NLPAN is specified using the parameters NMUSE, MHIN(I), NSOL(I), ISOL(U), FORCE, and SYMSTR. 
Features for reviewing a range of VIPASA eigensolutions are available in order to help select which 
modes to use. Initially, PASCO is called and instructed to generate the primary buckling solution for each 
longitudinal halfwave number from 1 to MSRCH. The eigenvalues are printed to ’nlpan.out’ along with 
an indication of whether or not the mode shape is symmetric (based on input parameters SYMSTR and 
NODSYM(l-2) in ’nlpan.in’.) To inspect these results before further execution, set IHALT=1 in ’nlpan.in’. 
In order to investigate the characteristics of secondary buckling eigensolutions for selected halfwave 
numbers, set IHALT=2 in ’nlpan.in’. NMUSE and MHIN(I) are used to specify which halfwave numbers 
are explored, and for each halfwave number MHIN(I), the number of solutions generated is equal to twice 
the maximum index number specified in ISOL(I,NSOL(I)). Both eigenvalues and symmetry indicators are 
printed to ’nlpan.out’ for all modes generated, and program execution is halted. 

Users may want to set up strategies for automatically selecting mode sets based on particular criteria. 
The appropriate place for this is subroutine PRIMAR in file ’primar.f’. A strategy implemented by 
Christine Perry for a certain class of problems can be found there and perhaps modified to suit the user. 

Restart procedure. The bulk of the computational effort of NLPAN (in the absence of extensive post- 
processing) goes into the stages which culminate in the computation of the coefficients of the nonlinear 
algebraic equations (stages (1) through (4) listed above). The user may choose to have the program store 
the results of all calculations made up to this point. This allows an unlimited number of subsequent runs 
to be made (using the stored information) in which the user may vary the shape imperfection, the load step 
sizes, the asynchronous load strategy, and other parameters affecting the nonlinear solution strategies. The 
operating mode is selected by answering the following prompt which appears on the screen at the start 
of the program execution: 

ENTER: 1 - NEW START WITH NO SAVE, 2 - NEW START WITH SAVE, 3 - RESTART 
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The RESTART mode does not restart a sequence of load/response calculations, but rather begins the 
computation of one or more new load-response sequences. In RESTART mode, the program reads the 
saved results, followed by input file ’nlpan.in2\ Output to file ’nlpan.out’ begins directly with a report 
of the final nonlinear analysis. 

If it is desired to suppress the screen prompts, then replace the statement 199=6 with 199=99 in 
subroutine NLP AN in source code file ’nlpan.f and recompile the program. The restart feature will not 
be used. Screen-directed output will be written to file ’term.out’. 

Special output. During each NLP AN run, output file ’nlpan.out’ is written which contains a summary 
of input information, properties of the VIPASA buckling modes, and other basic information about the 
model and the analysis to be performed. Basic output for the equilibrium solution at each load step 
includes normalized global load measures, modal amplitudes, and a reference deflection value. Output file 
’nlpan.outr’ contains a report, at each load step, of the values of various parameters used in the nonlinear 
solution strategies. (Definitions of the parameters listed in ’nlpan.outr’ can be found in comment lines in 
the beginning of source-code file ’nlanl.f .) Additional output is available for various stages of the analysis 
using print flags ILNPRT, ICOORD, IEFPRT, IQUPRT, INGPRT, and ICOPRT in file ’nlpan.inl’, and 
post-processing flags IPDISP, IPSTRN, IPSS, IPQY, IPSRES, and IPPROF in ’nlpan.in2’. The details of 
using these flags are described in the section "Line-By-Line Input Instructions." 

NLP AN is also capable of creating configuration, buckling mode, and displacement output files which 
can be viewed graphically using the commercially available PATRAN computer software by PDA 
Engineering of Costa Mesa, California. To create these output files, ICOORP must be set to 2 in file 
’nlpan.inl’. To create the files containing buckling mode-shapes, set IEFPRT=3 in file ’nlpan.inl’. To 
create displacement output files, set IPPROF=l in file ’nlpan.in2\ Additional details regarding the 
PATRAN-readable files are included in the descriptions of parameters ICOORP, IEFPRT, NXINT, NEY, 
and IPPROF, found in the section "Line-By-Line Input Instructions." An example of a PATRAN-generated 
image of displacement results from an NLPAN analysis is presented in Figure E2 for the case of a unit- 
stiffener-cell representation of a T-stiffened panel. 

Convergence Considerations 

The cross section of a configuration is discretized in order to perform numerical calculations for some 
steps of the NLPAN analysis. This discretization is determined by the parameter NINTN or parameters 
NINT(I) in ’nlpan.inl’. Increasing the fineness of the discretization improves the accuracy of results, but 
increases computer memory requirements and increases program execution time. The convergence of 
results with respect to the fineness of the discretization should be considered when using the program. 
Convergence is slowest when the transverse in-plane load component NY is non-zero. 

The accuracy of the NLPAN analysis is highly dependent on the suitability of the VIPASA buckling 
mode set used to represent displacements in the nonlinear analysis. NLPAN program features are provided 
to help the user screen the buckling modes for suitability, but the question of which modes to include is 
a difficult one. Mode-selection strategies for some cases are discussed in Section 5.5.2. 

Possible Problems 

VIPAS A/PASCO returns buckling mode shapes in terms of the associated generalized displacements 
of the node lines of the configuration. The mode shapes are normalized by setting the largest node-line 
rotation to 0. 1 . For some buckling modes (column type modes) all node-line rotations are very small, or, 
in theory, zero. This means that the associated displacements would be infinite in comparison. Experience 
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Figure E2 . Example of the visualization of NLPAN displacement results using PA TRAN, 
for a unit-stiffener-cell representation of a T-stiffened panel. 
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has shown that in practice, there is generally enough numerical error so that displacements are finite, 
though they may have amplitudes with an order of magnitude lxlO 6 * * * 10 . Mode shapes which have been 
plotted for such cases look physically reasonable. 

The solution strategies used to solve the nonlinear algebraic equations require the specification of 
several parameters to determine step sizes and cutoff values used in strategy decisions. In regions of severe 
solution-path curvature or compound solution branching, the solution strategies sometimes misdiagnose 
the situation, and issue a confusing series of messages. These problems can generally be eliminated by 
adjusting the values of DSNOM, DSFACT, DSMIN, DELQ, or EVCUT in ’nlpan.in2’. As far as what 
values to use, some trial and error may be required to obtain meaningful behavior. 

Avoid using extremely small non-zero values ( < 0.001 ) of modal imperfection amplitudes Q0(I) 
(specified in ’nlpan.in2’). Such values tend to create problems of the sort described in the previous 
paragraph. 

Line-By-Line Input Instructions 

In this section, detailed line-by-line instructions are given for the two input files ’nlpan.inl’ and 
’nlpan.in2\ The detailed instructions for input file ’pasco.in’ are covered in [El]. See the section "General 
Aspects of Specifying Geometry" in this appendix for special requirements which ’pasco.in’ must satisfy. 

General iastructions. The following general comments apply to the line-by-line input instructions, or 
to the input files. 

1) The contents of the input files must match, line per line, the data sequences described in this section. 

2) All input is read by FORTRAN as unformatted input. All input variables must be present, with at least 
one space separating data entries. Integer, floating-point, character, and logical input values must be 
written using conventional FORTRAN formats. Characters following the last required input value on a 
line are ignored. 

3) After all the required input values on a line have been read, the rest of the line is ignored, so that 
remaining space can be used for writing comments. 

4) Any special output requested using input parameters is written to the standard output file ’nlpan.out’, 
except where noted otherwise in the instructions. 

5) All CHARACTER input values must be enclosed by apostrophes except for the first record (TITLE). 
Example: SYMSTR may be ’Y’ or ’N’ . 

6) The variable type for each input item is indicated under the variable name, using the following 

abbreviations: 

I - Integer 

R - Real 

C*n - Character variable of length n 
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File ’nlpan.inl’ input lines. (This file is used to specify geometry, and modelling and output options.) 


TITLE 

C*76 

TITLE 


ICOORP 

I 

ICOORP 


ILNPRT 


IEFPRT 


A descriptive title used for labeling of the output files. 


ILNPRT IEFPRT IQUPRT 

I I I 

Controls printing of the global y-z coordinates of the model at the discretization 
points 

0 No print 

1 Print the Y-Z coordinates of the cross section to file ’coords.out’ 

2 Create file ’nlpan.neu’ which contains a PATRAN-readable "neutral file" describing 
the geometric configuration. This file contains finite-element node and element 
definition records, allowing the configuration to be viewed graphically using the 
PATRAN computer software. Discretization of the configuration is determined by 
input parameters NX ENT and NEY(IP) (IP=1,2,...,NPLATE). The file-creation options 
for writing PATRAN-compatible neutral results files for buckling mode shapes and 
displacements are enabled (see IEFPRT in file ’nlpan.inl’ and IPPROF in file 
’nlpan.in2’). 

Controls the operation of the linear solution module. See the discussion in the section 
"NLP AN Modes of Operation." 

0 Use the global unit loads NX(1) and NY(1) specified in ’pasco.in’, and force ILSET 
to 1 

2 Pre-processor mode. Based on the ILSET input value, select the two generalized in- 
plane load components from the values NXGL, NYGL, EPXL and EPYLG 
(described below); write corresponding PASCO input values NX(1) and NY(1) (and, 
if COMPAT=’Y\ any non-zero PASCO input factors FNY(I)) to the file ’nlpan.out’, 
for later insertion into ’pasco.in’ using a text editor. Program execution is halted. 

Controls the printing of buckling mode shapes, measured with respect to local plate- 
strip coordinate directions except where noted otherwise. 

0 No print. 

1 Print displacements. 

2 Print displacements and all computed derivatives. 

3 If ICOORIV2, then print displacements to file ’modes.out’ and halt execution. 

If ICOORP=2, then print mode shapes to PATRAN-readable files ’nlpanMMM.eig’ 
(MMM=1,2,...) where MMM is the index number of the buckling mode. These are 
"neutral results files" (in PATRAN nomenclature) containing buckling-mode 
displacements, and are to be used in conjunction with the configuration neutral file 
’nlpan.neu’ for graphical viewing of the buckling mode shapes. The displacements 
are measured with respect to global coordinate directions. Program execution is 
halted. 
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IQUPRT Controls printing of second-order displacement fields at discretization points, 

measured with respect to local plate-strip coordinate directions. 

0 No print. 

1 Print displacements. 

2 Print displacements and all computed derivatives. 

3 Print displacements to file ’uij.out’ and halt execution. 

INGPRT ICOPRT IPRINT IPBUG NXINT 

I I I I I 

INGPRT Controls the printing of the coefficients for the expression of global loads NXG and 

NYG (computed in subroutine ACALC). 

0 No print 

1 Print values 

ICOPRT Controls the printing of the coefficients of the nonlinear algebraic equations 

(computed in subroutine CCALC). 

0 No print 

1 Print coefficients 

2 Print coefficients and primitive coefficients 

IPRINT Development use only. Set to 0 

IPBUG Development use only. Set to 0 

NXINT Required only if ICOORP=2 (creating PATRAN-readable output files). NXINT is the 

number of finite elements along a length-wise cut of any plate strip in the structure, 
used in producing PATRAN-readable files ’nlpan.neu’, ’nlpanMMM.eig’, and 
’nlpanNNN.dis’. This is used for graphical visualization purposes only, and does not 
affect the analytical results. 


USER IHALT MSRCH 

C*5 I I 

USER Flag used by program developers. Set to ’STOLL’. 

IHALT Used to control the generation and display of the VIPASA buckling eigensolutions, 

and the subsequent use of the buckling mode shapes in the nonlinear analysis. See 
the discussion in the section "NLP AN Modes of Operation." 

0 Proceed with the NLP AN analysis using buckling modes as determined by the input 
values NMUSE, FORCE, MHIN(I), NSOL(I), ISOL(I,J) 

1 Compute the VIPASA buckling eigenvalues for longitudinal-halfwave numbers 1 
through MSRCH; print the eigenvalues along with an indication of the symmetry of 
the buckling mode. Halt execution. 

2 For values of the halfwave number specified by MHIN(I) (1= 1 ,2 NMUSE) compute 

and print several VIPASA eigenvalues, including an indication of the symmetry of 
each buckling mode. The number of modes computed for MHIN(I) is 
2*ISOL(I,(NSOL(I)). Halt execution. 
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MSRCH 


This parameter determines the range of longitudinal halfwave numbers investigated 
when IHALT=1 (see above). MSRCH must be in the range 0 to 20 
0 Investigate halfwave numbers 1 through 20 
>0 Investigate halfwave numbers 1 through MSRCH 


SYMSTR NODSYM(l) NODSYM(2) 

C*1 I I 

SYMSTR ’Y’/’N’ Indicates whether or not the structural cross section is symmetric. 

NODSYM(J) Index numbers of two node lines for which global displacements W can be compared 
to detect symmetry of a buckling mode. Used only if SYMSTR=’Y’ 


NMUSE FORCE 

I C*1 

These parameters, along with parameters NSOL, MHIN, and ISOL in following records, 
determine which buckling modes (computed by VIPASA) are used in the nonlinear 
analysis. 

NMUSE Number of unique values of the longitudinal halfwave-number for which VIPASA 

mode shapes will be generated for use. 

FORCE Affects the selection of VIPASA buckling modes for the longitudinal halfwave- 

numbers specified with MHIN(I). In the following description, designate buckling 
solutions by the indices (M,J). Solution (M,J) is the J’th solution in the infinite 
sequence of solutions (for halfwave number M) which are ordered in terms of 
eigenvalue. 

’Y’ For longitudinal halfwave number MHIN(I), force the use of modes 
(MHIN(I),ISOL(I,J)), J= 1 ,2 NSOL®. 

’N’ Input values ISOL are ignored. Use mode (MHIN(I),1). If SYMSTR=’Y’ and 
NS0L(I)>1 then add the next NSOL(I)-l modes (MHIN(I),N) which match mode 
(MNIH(I),1) with respect to symmetry or lack thereof. If SYMSTR=’N’ use modes 
(MHIN(I),J), J=l,2,...JMSOL(I). 


NSOL(l) NSOL(2) ... NSOL(NMUSE) 

I I I 

NSOL(I) Number of VIPASA mode shapes to be incorporated which have the longitudinal 

halfwave-number MHIN® 


MHIN(l) MHIN(2) ... MHIN(NMUSE) 

I I I 

MHIN® Longitudinal halfwave-number of a VIPASA mode shape to be used in the NLP AN 

analysis 
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IS0L(1,1) IS0L(1,2) ... IS0L(1,NS0L(1)) 

IS0L(2,1) ISOL(2,2) ... ISOL(2,NSOL(2)) 

ISOL(NMUSE, 1 ) IS0L(NMUSE,2) ... ISOL(NMUSE,NSOL(NMUSE)) 

I I I 

ISOL(IJ) Indicates which VIPASA buckling modes to use from the sequence of modes having 
longitudinal halfwave-number MHIN(I), where modes are ordered according to 
eigenvalue. All input values must be present, but they are only used if FORCE=’Y’. 


NPLATE NNODE 

I I 


NPLATE Number of plate strips in the model 
NNODE Number of node lines in the model 


NPOFFS 

I 

NPOFFS Number of plate strips which have non-zero offsets between one (or both) side edge(s) 
and the corresponding node line(s). 


Conditional - include only if NPOFFS >= 1 : 

IP 1ECVY(IP,1) IECVZ(IP,1) IECVY(IP,2) IECVZ(IP.2) (1=1) 
IP IECVY(IP.l) IECVZ(IP,1) IECVY(IP,2) IECVZ(IP,2) (1=2) 


IP IECVY(IP.l) IECVZ(IP,1) IECVY(IP,2) 

II I I 


IECVZ(IP,2) 

I 


(I=NPOFFS) 


IP Index number of a plate strip with non-zero node-line offsets. 

IECVY(IP,IE),IECVZ(IP,IE): Integer values used to determine the offset component values. This is 
done the same way as in PASCO. The absolute value gives the index number I of a 
thickness value H(I) specified in ’pasco.in’, and the sign specifies the direction of the 
offset component. IE is the edge number, IE=1 for Y=0, IE=2 for Y=B(IP). 
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MUA(IP) SKVEC(IP) PRVEC(IP) 

MUA(IP) SKVEC(IP) PRVEC(IP) 


TUA(IP) (IP=1) 
TUA(IP) (IP=2) 


MUA(IP) 

R 


SKVEC(IP) PRVEC(IP) TUA(IP) (IP=NPLATE) 

I I R 


MUA(IP) Rotational orientation angle, in degrees, of plate strip # IP in Y-Z plane. Same 

convention as PASCO. See Figure El(b). 


SKVEC(IP) Indicator for presence of an initially flat skin which is continuous between the 
boundary node-lines. 

1 Plate strip # IP is part of the skin. For the case of overlapping plate strips, such as 
a skin/stiffener-flange region, only one of the overlapping strips should be assigned 
SKVEC(IP)=1 

0 Plate strip # IP lies off of the skin. 


PRVEC(IP) 


0 

1 

-1 


Used for indicating pressure load application and direction. Note: For the case of 
overlapping plate strips, such as a skin/stiffener-flange region, only one of the 
overlapping strips should be assigned a non-zero value PRVEC(IP) 

No pressure load applied to plate strip 

Pressure load acts on plate strip in local +Z direction. 

Pressure load acts on plate strip in local -Z direction. 


TUA(IP) Unit value of temperature (the difference from a reference temperature) to be applied 

uniformly to plate strip # IP. See the description of parameter HEATA in input file 
’nlpan.in2’. 


NOD(IP,l) NOD(IP,2) (IP=1) 

NOD(IP,l) NOD(IP,2) (IP=2) 

NOD(IP,l) NOD(IP,2) (IP=NPLATE) 

I I 

NOD(IP.IE) Node-line index number of edge # IE of plate strip # IP. Defines the connectivity of 
the model by specifying the two node lines to which each plate strip attaches. Node- 
line numbers and plate strip numbers must match those used to specify the geometry 
in ’pasco.in’. Plate-strip numbers are sequential, and correspond to initial PASCO 
plate strip numbers (i.e. before the application of HCARD conversions). IE is the 
edge number, IE=1 for Y=0, IE=2 for Y=B(IP). 


BNODE(l) BNODE(2) 

I I 

BNODE(IB) Index numbers of two "boundary" node lines at which boundary conditions are applied. 

IB=1 nominally corresponds to the global Y=0 boundary, IB=2 nominally corresponds 
to the Y=B(Global) boundary. The boundary node lines need not span the width of the 
configuration, but if they do not, then global load NY may not be computed correctly. 
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BCVEC(l.IB) BCVEC(2,IB) BCVEC(3,IB) BCVEC(4,IB) (IB=1) 

BCVEC(l.IB) BCVEC(2,IB) BCVEC(3,IB) BCVEC(4,IB) (IB=2) 

I I I I 

BCVEC(IC.IB) Boundary-condition indicator for component # IC of boundary node-line # IB. The 
following table lists options and components: 

IC Comp. BCVECUC.IB) 

1 U, Fx 1 or 2 

2 V, Fy 1 or 2 or 3 

3 W, Fz 1 or 2 

4 dW/dy, M 1 or 2 

Note that the four components are specified in the opposite order compared to 
PASCO. 

BCVEC = 1 Control the generalized displacement Set to zero (except for IC=2, where V may be 
non-zero). 

BCVEC = 2 Control the generalized force. Set to zero (except for IC=2 for which Fy (=NY) may 
be non-zero). 

BCVEC = 3 IC=2 only. Keep the edge straight w.r.t. global V, but control the mean global NY 
load. 


NCU NCW 

I I 

NCU Number of axial displacement constraints, each end. Axial displacement U is constrained 
to follow the effective end-shortening. This is imposed only on modes with odd 
longitudinal halfwave numbers. These constraints may only be used with CONTRL=’D\ 

NCW Number of slope (dW/dX) constraints, each end. The slope is measured with respect to 
the local coordinate directions of specified plate strips. 


Conditional - include only if NCU > 1 : 

IPCU(l) YCU(l) 

IPCU(2) YCU(2) 

IPCU(NCU) YCU(NCU) 

I R 

IPCU(I) Plate strip on which the axial displacement constraint is imposed. 

YCU(I) Nominal location Y on plate strip # IPCU(I) where the constraint is applied. The actual 

point of application will be at the closest available discretization point. See the description 
of parameters NINTN and NINT(IP). 
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Conditional - include only if NCW £ 1 : 

IPCW(l) YCW(l) 

IPCW(2) YCW(2) 

1PCW(NCW) YCW(NCW) 

I R 

IPCW(I) Plate strip on which the slope displacement constraint is imposed 

YCW(I) Nominal location Y on plate strip # IPCW(I) where constraint is applied. The actual point 

of application will be at the closest available discretization point. See the description of 
parameters NINTN and NINT(IP). 


NSPW 

I 

Number of rotational springs (at each end) restraining dU/dY (Y - Local plate-strip axis). 

Number of rotational springs (at each end) restraining dW/dX. (W - Local plate-strip 
deflection). 


Conditional - include only if NSPU > 1 : 

IPSPU(l) YSPU(l) KSPU(l) 

IPSPU(2) YSPU(2) KSPU(2) 

IPSPU(NSPU) YSPU(NSPU) KSPU(NSPU) 

I R R 

IPSPU(I) Plate strip on which rotational spring acts (each end) 

YSPU(I) Nominal location Y on plate strip # IPSPU(I) where the spring acts. The actual location 
will be at the closest available discretization point. See the discussion of discretization 
using parameters NINTN and NINT(IP). 

KSPU(I) Rotational-spring constant (Moment/Radian) 


Conditional - include only if NSPW > 1 : 

IPSPW(l) YSPW(l) KSPW(l) 

IPSPW(2) YSPW(2) KSPW(2) 

IPSP W (N SP W) YSPW(NSPW) KSPW(NSPW) 

I R R 

IPSPW(I) Plate strip on which rotational spring acts (each end) 


NSPU 

I 

NSPU 

NSPW 
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YSPW(I) 


Nominal location Y on plate strip # IPSPW(I) where the spring acts. The actual location 
will be at the closest available discretization point. See the discussion of discretization 
using parameters NINTN and NINT(IP). 

KSPW(I) Rotational-spring constant (Moment/Radian) 


Comment line - not used for data. 

**** DISCRETIZATION AND REFERENCE VALUES **** 


NINTN 

I 

NINTN Guides the discretization of each plate strip in the local Y-direction. Used for the 

finite-difference computation of the second-order displacement fields and for 
numerical integration. 

0 Specify the discretization of each plate strip individually using NINT(IP) 

N N is an even integer >4. The widest plate strip in the structure will have N+l 
discretization points evenly spaced over N intervals. The remaining plate strips are 
discretized such that interval widths on all plate strip are approximately equal, except 
that no plate strip will have less than 4 intervals. 


Conditional - include only if NINTN = 0 : 

NINT(1) NINT(2) ... NINT(NPLATE) 

I I I 

NINT(IP) Even integer > 4 giving the number of discretization intervals in the Y-direction on plate 
strip # IP. Note: Plate strips with matching indices ILBWAL(IP) must have matching 
numbers NINT(IP). 


Conditional - include only if ICOORP=2 (creating PATRAN-readable output files): 

NEY(l) NEY(2) ... NEY(NPLATE) 

I I I 

NEY(IP) Even integer > 2 giving the number of discretization intervals in the Y-direction on plate 
strip # IP to be used for graphical visualization with PATRAN. Generally, NINT(IP) is 
so large that if all available intervals are used in producing graphical images in PATRAN, 
the images are too cluttered. NEY(IP) must be selected so that NINT(IP) is an integer 
multiple of NEY(IP), i.e. NINT(IP)=N*NEY(IP), where N is a positive integer. 
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IPHREF IPWDET IPDFL IYDFL 

I I I I 

IPHREF Controls the reference thickness used in normalizing modal amplitudes and 

displacements. The reference thickness HREF used by NLP AN is listed in file 
’nlpan.out’. 

0 HREF is the thickness for plate strip # 1 
IP HREF is the thickness for plate strip # IP 

IPWDET Determines the plate strip where the characteristic amplitude of a buckling mode 

(used for normalizing the modal amplitudes) is measured 
0 The characteristic amplitude is the maximum displacement in the structure 
IP The characteristic value is the maximum displacement on plate strip # IP 

IPDFL IPDFL and IYDFL control the transverse location where a characteristic displacement 

is computed and reported during the nonlinear analysis. The displacement 
(normalized by HREF) is printed in ’nlpan.out’ for each equilibrium solution, labeled 
as WCN. This feature serves only as a user convenience. 

IP Deflection is measured on plate strip # IP 

IYDFL Determines the location on plate strip # IPDFL where the characteristic deflection is 

computed. Where B is the width of plate strip # IPDFL: 

-1 Compute deflection at Y=0, X=L/2, plate strip # IPDFL 

0 Compute deflection at Y=B/2, X=L/2, plate strip # IPDFL 

1 Compute deflection at Y=B, X=L/2, plate strip # IPDFL 


Comment line - not used for data. 

**** LOADING AND MODELLING OPTIONS **** 


CONTRL COMPAT 

C*1 C*1 

CONTRL Specifies the type of control used for the generalized in-plane loads. 

Note: for simple rectangular plates, both types of control yield the same load/end- 
shortening relationship, but for more complex configurations, a discrepancy arises 
between the two methods. A correction is under development to fix this 
discrepancy, but currently CONTRIV’D’ is believed to provide greater accuracy. 

’D’ Displacement control 

’L’ Load control 

COMPAT Specifies whether displacement compatibility is enforced in the linear solution for 

multiply connected cross sections. See the discussion in the section "NLP AN 
Modes of Operation." 

’Y’ Displacement compatibility is enforced. In general, NLP AN must be run first in 
a pre-processor mode (ILNPRT=2) to compute values NX(1), NY(1), FNY(IP) to 
put in file ’pasco.in’. These are printed in file ’nlpan.out’. 

’N’ Displacement compatibility is either not enforced, or is not applicable. 
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ILSET NXGL NYGL EPXL EPYLG 

I R R R R 


Note: This input set must be present, but values are used only if ILNPRT=2 (pre-processor 

mode). If ILNPRT=0 then ILSET=1 is assumed, and NXGL and NYGL (NX(1), NY(1» 
are taken from file 'pasco. in’. 


ILSET 


Determines which two parameters define the unit in-plane load system. All loads and 
strains defined positive for tension/extension. Unit in-plane loading must generally 
be tensile, because buckling eigenvalues are assumed to be negative. 

1 NXGL, NYGL 

2 NXGL, EPYLG 

3 EPXL, EPYLG 

4 EPXL, NYGL 


NXGL 

NYGL 

EPXL 


Positive unit value for global X-normal axial load per unit width based on global 
width B. Equivalent to NX(1) in PASCO input. 

Unit value for global Y-normal axial load per unit length based on length L. 
Equivalent to NY(1) in PASCO input. 

Unit value for axial strain in X-direction 


EPYLG 


Unit value for mean Y-normal strain in skin (change in panel width per unit width) 


PRUNIT 

R 

PRUNIT Unit pressure load, force per unit area. See the description of parameter BET A A in 

input file ’nlpan.in2’. 


LOCGLO 

C*1 

LOCGLO Y / N Used to control modifications to the second-order displacement fields necessary 
to improve the solution accuracy, particularly in cases of local-global mode 
interaction. In general, set LOCGLO=’Y’ unless INDPLT=1 (simple plate 
analysis). 


C-l 



INDPLT 

I 

INDPLT 1 Configuration is a rectangular plate with a single reference plane (may be constructed 
from several linked plate strips if no node-line offsets are used). Implies that 
u i =v i =w ij = 0- This reduces the cost of the analysis by avoiding unneeded calculations 
and data storage. 

2 Complex configuration with multiple reference planes. Implies that Uj, v-, w ; are in 
general not zero. 1 IJ 


Comment line - not used for data. 

**** LOAD, WIDTH, WALL INDEXING **** 


ILBWAL(l) ILBWAL(2) ... ILBWAL(NPLATE) 

I I I 


ILBWAL(IP) See the discussion in the section "General Aspects of Specifying Geometry." 

ILBWAL(IP) is an integer in the interval 1 to NLBWAL, where NLBWAL is the 
number of unique combination of: 

i) wall properties, 

ii) width B, and 

iii) unit in-plane loads, 

among all plate strips in the structure. ILBWAL(IP) specifies which set of 
characteristics plate strip # IP has. These indices must match the values assigned by 
PASCO. NLBWAL is passed to NLP AN via file ’modes.in’ created by PASCO. 


NKWALL 

I 

NKWALL Number of different laminates defined in input file ’pasco. in’. 
— End of file ’nlpan.inl’ — 
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File ’nlpan.in2’ input lines. (This file contains the control parameters for the nonlinear analysis, 
the control parameters for post-processing, and the modal-imperfection amplitudes. File ’nlpan.in2’ is read 
every time NLP AN is run, including runs in the RESTART mode.) 


NIMP IZCIP SWITCH 

I I Logical 


NIMP 

IZCIP 


Number of sets of modal imperfection amplitudes for use in an equal number of 
consecutive nonlinear analysis run. NIMP lies in the range 1 to 30 

Controls the zeroing of coefficients CIP(I) (in the nonlinear algebraic equations) 
which, when non-zero, cause the initial response to in-plane loading to be nonlinear. 
Each coefficient CIP(I) corresponds to a particular VIPASA buckling mode 
incorporated in the NLPAN analysis. CIP(I) may be non-zero either for physical 
reasons, indicating that the initial response is truly nonlinear, or may be non-zero due 
to approximations used in the method of analysis. One application is to set IZCIP=1 
to eliminate an initial bowing response characteristic of bowing imperfections. 

0 Do not zero the coefficients CIP(I). 

-1 Set CIP(I) to zero for all modes. 

M Set CIP(I) to zero if NWAVEA(I)<= M, where NWAVEA(I) is the longitudinal 
halfwave number for the VIPASA buckling mode 


SWITCH Flag used in an investigation of compound solution branching. Set to .FALSE.. (If 

set to .TRUE., then at points of near-simultaneous solution branching in two modes, 
the initial branch-switching step will be controlled by a specified increment to the 
non-critical mode rather than by a specified increment to the critical mode.) 


Q0(l) 

Q0(i) 

Q0(2) ... 
Q0<2) ... 

QO(NEIG) 

QO(NEIG) 

(IIM=1) 

(IIM=2) 

Q0(l) 

R 

Q0(2) ... 
R 

QO(NEIG) 

R 

(IIM=NIMP) 


Q0(I) Modal imperfection amplitude for VIPASA buckling mode number I, where I varies from 
1 to NEIG, where NEIG is the number of modes incorporated in the analysis. Modal 
amplitudes are normalized based on a maximum deflection equal to the reference 
thickness value HREF (listed in ’nlpan.out’). Characteristics of the modes including 
critical load level, halfwave number, symmetry, and total number NEIG are listed in the 
output file ’nlpan.out’. 
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QF(NEIG) 

R 


QF(1) QF(2) ... 

R R 

QF(I) Modal amplitudes used to specify forced end displacements or end rotations. QF input 
must be present, but is used only when displacement constraints are used (NCU > 2 or 
NCW > 1). These values differ from imperfection amplitudes QO(I) in that they are stress- 
producing. Amplitudes QF are significant only in terms of the associated generalized 
displacements at the longitudinal ends. They are intended for use in simulating a forced 
initial rotation of the panel ends. 


NRANGE 

I 


NRANGE Number of load ranges specified. Intended for use in modelling asynchronous application 
of multiple load types. Within each load range, the in-plane, pressure, and thermal load 
parameters all vary linearly with a single generalized load parameter. All generalized loads 
are continuous from one range to the next, but the load rates are discontinuous. NRANGE 
must be in the range 1 to 5. 


LAMDAA(IR) BETAA(IR) HEATA(IR) (IR=1) 

LAMDAA(IR) BETAA(IR) HEATA(IR) (IR=2) 


LAMDAA(IR) BETAA(IR) 
R R 


HEATA(IR) (1R=NRANGE) 
R 


These are the target values for the three load parameters, used to define the NRANGE 
load ranges. All three parameters are assumed to be zero initially. 


LAMDAA In-plane load parameter. This is a normalized value, where the reference value REFLAM 
is listed in output file ’nlpan.out’. LAMDAA=1.0 corresponds to the critical buckling load 
for the case of pure in-plane loading. 


BETA A Load multiplier for transverse pressure. The applied pressure P is given by 
P=BETAA*PRUN1T, where PRUNIT is the unit pressure specified in file ’nlpan.inl’. 

HEATA Load parameter for thermal loading. In a given plate strip, the thermal load T is given by 
T=HEATA*TUA(IP), where TUA(IP) is the unit temperature specified in file ’nlpan.inl \ 


91 



TOL ITMAX UU VV 

R I R R 

These arc input parameters for subroutine ZROOT, which uses Barstow’s method to find the roots of 
a polynomial with real coefficients. The polynomial equation arises when applying Thurston’s method 
for negotiating solution branching. The default values provided generally work well, and need to be 
changed only if a warning message is issued. 


TOL 

0 . 

Tolerance to measure convergence. Default = l.E-12 
Use default value. 

ITMAX 

0 

Maximum number of iterations allowed. Default = 20 
Use default value 

UU, VV 


Starting values used in iterative solution procedure. The starting values are 
automatically varied in NLP AN if convergence problems are indicated. Defaults: 
UU=-1., VV=2. 


0 . 

Use default value. 


ILSTM CCRIT MAXIT 

I R I 


ILSTM Positive integer giving the maximum permitted number of load steps. 

CCRIT Tolerance to measure convergence of the iterative procedures for solving the nonlinear 
equilibrium equations. Suggested value: l.E-5 

MAXIT Maximum number of iterations allowed in the iterative solutions procedures. With the 
methods used, convergence is usually rapid. Suggested value: 6 


DSNOM DSFACT DSMEN DELQ EVCUT DSQUIT 

R R R R R R 

These parameters control the solution step sizes, and specify cutoff and decision values in the nonlinear 
solution strategies. DS refers to the arc-length increment used from one solution to the next in a Riks- 
type approach of arc-length control. The modal amplitudes, load parameters, and Lagrange multipliers 
(when needed) are all normalized to take on values with order of magnitude unity, so DS should be 
interpreted accordingly. For all six parameters, an input value 0.0 forces the use of the default value. 

DSNOM Nominal value used for DS. DS will never exceed this value. Default: 0.3 

DSFACT Factor for reducing DS compared to DSCRIT, where the latter is the arc-length increment 
predicted (based on extrapolation) to lead to a critical stability point. If 
(DSFACT*DSCRIT)<DSNOM then DS=DSFACT*DSCRIT. Default: 0.6 

DSMIN Minimum value permitted for DS before a solution at the critical point solution is 
attempted. Default: 0.002 
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DELQ Step size used for initiating a solution along a new path at a branch point. Approximately 
equivalent to an increment in a modal amplitude. Default: 0.005 

EVCUT If a solution branch point is detected, EVCUT is used to select between a one-mode and 
two-mode solution branching analysis. If the first two eigenvalues of the tangent stiffness 
matrix are separated by more than EVCUT, a single-mode branching analysis is 
performed, otherwise a two-mode branching analysis is performed. Setting EVCUT small 
tends to force a single-mode branching analysis, in which case DELQ must also be small 
to avoid stepping past secondary branch points. Default: 0.01 

DSQUIT Minimum permitted value for DS. If DS is less than DSQUIT and solution convergence 
is not achieved, execution is halted. Default: 0.24*DSMIN 


IFREQ 


IFREQ >1 Postprocessing requests are processed at every IFREQ ’th load step. 
0 No post-processing 


(Note: This and all following input records are required only if IFREQ>0) 


IPDISP 

IPSTRN 

IPSS 

IPSRES 

IPQY 

IPPROF 

I 

I 

I 

I 

I 

I 


Post-processing control parameters. (Results are computed at model locations which 
are specified in following input records.) Each parameter may be specified as 0 or 
1, with the following effect: 

1 Output is requested. 

0 No output is requested 

IPDISP Print displacement values to file ’displ.out’ 

IPSTRN Print mid-surface strains and curvatures to file ’strain.out’ 

IPSS Print surface strains EPSx and EPSy to file ’strain.out’ 

IPSRES Print stress resultant values to file ’sres.out’ 

IPQY Print force resultant QY which acts in the Z-direction on the Y-normal face in a plate 

strip, to file ’sres.out’ 

IPPROF If ICOORP*2: print the displacement profile (referred ro local plate-strip coordinate 

directions) at a constant X-station to file ’displ.out’. 

If ICOORP=2: print a PATRAN-readable displacement-results file ’nlpanNNN.dis’, 
where NNN is the corresponding load step number. The displacements are computed 
at locations determined by parameters NXINT and NEY(IP) of input file ’nlpaainl ’. 
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NINTX NXLOC 

I I 


NINTX >1 If NINTX>1 then NINTX is the number of uniform intervals in the X-direction 
determining where postprocessing results are computed. NINTX may be up to 100 
for displacements, and up to 50 for strains and stress resultants. If NINTX>1 then 
NXLOC is forced to NINTX+1. 

0 Use NXLOC and XLOCN input to specify X locations. 


NXLOC >1 Number of specific X values used for post-processing. Ignored if NINTX>1. 


Conditional - include only if NINTX=0 and NXLOC>0: 

XLOCN(IX) (IX=1) 

XLOCN(IX) (IX=2) 

XLOCN(IX) (IX=NXLOC) 

R 

XLOCN(IX) Value of X, as a fraction of length L, where results are to be computed. 0.0 < 
XLOCN(IX) < 1.0 


NYLOC 

I 

NYLOC Number of cross-sectional stations where results are computed. Maximum allowable: 100 


IP(IY) 

ID(IY) 

(IY=1) 

IP(IY) 

ID(IY) 

(IY=2) 

IP(IY) 

ID(IY) 

(IY=NYLOC) 


I I 


IP(IY) Index number of the plate strip where results are to be computed. 

ID(IY) Index number of the discretization point in the Y-domain of plate strip # IP(IY) where 

results are to be computed. Index ID(IY) lies in the range 1 to (MNT(IP)+1). Values 

NINT(IP) (IP=1,2 NPLATE) are listed in the output file ’nlpan.out’, and can be 

specified directly for each plate strip in input file ’nlpan.inl’. 

— End of file ’input.in2’ — 
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Example Input Files 


Unit cell representation of a T-stiffened panel. The input files for a test case designated UNITA are 
listed below. UNITA refers to a unit-cell representation of a T-stiffened panel. The general configuration 
is that shown in Figure El(a). This model was used to analyze a panel which was tested experimentally. 
A comparison of analytical and experimental results is presented in [E2], The test panel had three evenly 
spaced T-stiffeners. The analytical model includes only a single stiffener with a half-bay of skin on each 
side, with symmetry conditions imposed at the side edges of the skin. Three different laminates were 
defined to model the skin, the stiffener flange, and the stiffener blade, respectively, and the mean values 
of lamina thickness measured for the test panel were used in the analysis. The thickness of the adhesive 
used to bond the stiffeners to the skin is represented in file ’pasco.in’ by T(10); the in-plane stiffness of 
the adhesive is neglected. Thickness T(ll) is used to specify the offset between the mid-surface of the 
stiffener flanges and the mid-surface of the skin. Thickness T(ll) is defined in ’pasco.in’ in terms of the 
other thicknesses using a constraint equation established by parameters AT and AC. The offset values are 
assigned using the first two HCARD records in file ’pasco.in’, and using the IECVZ parameters in file 
’nlpan.inl’. 

Input file ’nlpan.inl’ features axial displacement constraints at the longitudinal ends at node-line 
numbers 3 and 4 (see Figure El (a)) in order to simulate clamping of the panel ends. The constraint 
locations are specified in terms of locations on plate strips 3 and 7 using IPCU and YCU input. The 
following modes, specified by input parameters NMUSE, NSOL, MHIN, and ISOL, are used in the 
analysis: (1,1), (1,3), (3,2), (3,3), (5,2), (5,3), (7,2), and (7,4). All of these modes have symmetric cross- 
sections. The symmetry of the modes is indicated in the output file ’nlpan.out’, and is determined by 
comparing the deflections of node-line numbers 2 and 5 (see the NODSYM input in *nlpan.inl’). 

In input file ’input.in2’, an imperfection amplitude Q0(5)=0.02 is specified (Q(5) corresponds to mode 
(5,2), the designated critical local mode), and forced end rotation is imposed through input QF(1)=1.79. 
Mode 1 (mode (1,1)) is a bowing type mode. The longitudinal variation of displacements for mode 1 is 
given approximately by [HREF SIN(PI X / L)], where HREF=0.0435 in. (the skin thickness) and L=20 
in. (the panel length). Using these values, it can be shown that QF(1)=1.79 corresponds to an end rotation 
of 0.7 degrees, which is the value of end rotation imposed on each end of the test panel during one 
experiment. Additional discussion of the analysis and tests is found in [E2], The PATRAN-generated 
image of a deformed panel presented in Figure E2 corresponds to this NLP AN model. The deformation 
state of the panel was computed for an equilibrium solution at the theoretical elastic limit-load point; the 
displacements have been amplified for clarity. 

In the listings of files ’nlpan.inl’ and ’nlpan.in2’, comments are set apart from input data using 
exclamation points (!). The exclamation points are used only for visual purposes, and serve no other 
function in the input files. 
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File ’pasco.in’ for the UNIT A example. 


UNIT A - TEST PANEL A - UNIT-CELL REPRESENTATION 

&CONDAT 

&END 

&PANEL 

B=1 .43,.49,.49, 1 .43,.49,.49,.49 I 
EL=20., 

CONV1=50000„ 

FREQ=0„ 

IP=0, 

MAXJJJ=0, 

MAT(1,1)=1,1, 1,2,2, 2,3, 3,3, 

T(l)=.00544,.00544,.00544,.00563„ 00563, .00563, 
.00483,.00483,.00483,.0045, 1 .E30, 
AT(1,1)=4.,0.,0.,4.,0.,0.,0.,0.,0.,1.,-1., 

AC(1)=0„ 

THET( 1 )=0, 50,90,0, 1 5 ,90,0, 1 5,90,0, 

KWALL(1 , 1 )=2,-2, 1 ,3, 

KWALL(l,2)=5,-5,4,6, 

KWALL(1 ,3)=8,-8,7,9, 

IWALL(1)=1,1,1,1,2,2,3, 

FSTIFF=10., 

HCARD=6,-8,-5,0,-l 1.0.-1 1, 

6,-9,-6, 0,-11, 0,-11, 

4,-10,7,90,0, 

2 , 121 , 10 , 

13,11,1,-990,9000,2,-8,3,-9,-121,4,0,-990,9000, 

ICARD=5, 1 ,2, 1 ,-990,9000, 

5.2.4.2.4.8, 

3,3,4,10, 

5.4.5. 3.5. 9, 

3.5, 6,4, 

3.6, -990,9000, 

ICREP=1, 

NOBAY=l, 

NX(1)=500., 

&END 

&MATER 

E1(1)=17.02E6, E2(1)=1.64E6, E12(l)=. 80E6, ANU1(1)=.30, 
RHO( 1)= 1 .477E-4, ALFAl(l)=0.25E-6, ALFA2(l)=16.2E-6, 
ALLO W( 1 , 1 )=2 ..0 1 1 , -0.0 1 1 ,0.0 1 ,-0.0 1 3 ,0.0 1 55 , 

E1(2)=16.44E6, E2(2)=1.64E6, E12(2)=.77E6, ANU1(2)=.30, 
RHO(2)= 1 .477E-4 , ALFA1(2)=0„ ALFA2(2)=0„ 

ALLO W(1 ,2)=2,.0 1 1 ,-0.0 1 1 ,0.0 1 ,-0.0 1 3 ,0.0 1 55, 

E1(3)=19.17E6, E2(3)=1.64E6, E12(3)=.90E6, ANU1(3)=.30, 
RHO(3)= 1 .477E-4, ALFA1(3)=0„ ALFA2(3)=0„ 

ALLO W( 1 ,3)=2,.0 1 1 ,-0.0 1 1 ,0.0 1 ,-0.0 1 3,0.0 155, 

&END 

— End of file ’pasco.in’ — 
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File ’nlpan.ini’ for the UNIT A example. 


0 


UNIT A - TEST PANEL A - UNIT-CELL REPRESENTATION 
ICOORP JLNPRT,IEFPRT,IQUPRT 
INGPRT 4COPRT,IPRINT,IPBUG 
USER, IHALT, MSRCH 
SYMSTR, NODSYM(l), NODSYM(2) 

NMUSE, FORCE 
NSOLfl), 1=1, NMUSE 
!MHIN(I), 1=1, NMUSE 
ISOL(U), J=1JMS0L(I), 1=1 
1=2 
1=3 
1=4 

INPLATE NNODE 
NPOFFS 

!IP,IECVY(IP,1),IECVZ(IP,1),IECVY(IP,2),IECVZ(IP,2) IPOFFS 

IPOFFS=2 

MU(IP), SKVEC(IP), PRVEC(IP) TUA(IP), IP=1 

IP=2 
IP=3 
IP=4 
IP=5 
IP=6 
IP=7 

NOD(IP,IE), IP=1 
IP=2 
IP=3 
IP=4 
IP=5 
IP=6 
IP=7 

BNODE(IB), IB=1,2 
BCVEC(IC,IB), IC=1 to 4, IB=1 

IB=2 

!NCU, NCW 

IPCU(ICU), YCU(ICU), ICU=1 
ICU=2 

NSPU, NSPW 

**** DESCRETIZATION AND REFERENCE VALUES **** 

ININTN 

IPHREF, IPWDET, IPDFL, IYDFL 
AND MODELLING OPTIONS **** 

CONTRL, COMPAT 
ILSET,NXGL,NYGL,EPXL,EPYLG 
PRUNIT 
ILOCGLO 
INDPLT 


0000 
0000 
’STOLL’ 0 
’Y’ 2 5 
4 ’Y’ 

2222 
13 5 7 

1 3 
23 

23 

24 
76 

2 

50-11 0-11 
60-110 -11 
0 . 1 1 0 . 

0 . 1 1 0 . 

0 . 1 1 0 . 

0 . 1 1 0 . 

0.0 0 0 . 

0.0 0 0 . 

90. 0 0 0. 

1 2 
24 

4 5 

5 6 
24 
4 5 
3 4 
1 6 

2 3 2 1 

2 3 2 1 
20 

3 0. 

7 0. 

00 


12 

102 1 
**** LOADING 
’D’ *N’ 

1 500. 0. 0. 0. 
0.0 
’Y’ 

2 


1221334 
3 


**** LOAD, WIDTH, WALL INDEXING **** 


ILB W AL(IP),IP= 1 NPLATE 
NKWALL 


- End of file ’nlpan.inl’ -- 
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File ’nlpan.in2’for the UNITA example. 


1 0 .FALSE. INIMP, IZCIP, SWITCH 

0. 0.0.0.0.02 0.0.0. !Q0NST(IE,IIMP),IE=1,NEIG) 

1.79 0. 0. 0. 0. 0. 0. 0. ! QF(IE) ,IE= 1 ,NEIG 

1 INRANGE 

20. 0. 0. !LAMDAA(1), BETAA(l), HEATA(l) 

0. 0 0. 0. !TOL UMAX UU VV 

30 0.00001 6 IILSTM, CCRIT, MAXIT 


.25 .6 .001 .005 .05 0. !DSNOM,DSFACT,DSMIN,DELQ,EVCUT,DSQUIT 

0 IIFREQ 

0 1 0 0 0 0 IIPDISP, IPSTRN, IPSS, IPSRES, IPQY, IPPROF 

0 3 ! NINTX.NXLOC 

0.094 !XLOCN(IX), IX=1 

0.406 ! IX=2 

0.5 ! IX=3 

2 INYLOC 

1 1 !IP(IY), ID(IY), IY=1 

7 3 ! IY=12 

— End of file ’nlpan.in2’ — 


Unit cell representation of a hat-stiffened panel. Input files ’pasco.in’ and ’nlpan.inl ’ are listed below 
for a unit-cell representation of a hat-stiffened laminated composite panel. The hat stiffener includes 
mounting flanges. The cross-sectional representation is shown in Figure E3. The slope of the webs of the 
hat stiffener (plate strip numbers 7 and 9) relative to the skin is 65 degrees for the configuration described 
in the input files. The mid-surfaces of the stiffener flanges are offset from the mid-surface of the skin by 
a distance d. Formulas for computing the offset measures EY and EZ between the edges of plate strips 
7 and 9 and the mid-surface of the skin are included in Figure E3, along with sketches which show the 
geometric details. The cross-sectional proportions and laminate configurations used in the input files were 
chosen arbitrarily. 

There are six thickness values T(I) included in input file ’pasco.in’. The first three values are lamina 
thicknesses. T(4) is the offset between the reference surfaces of the skin and the stiffener flanges. T(5) 
and T(6) are the amplitudes of the eccentricity measures EY and EZ, respectively, for plate strip numbers 
7 and 9, which are defined as illustrated in Figure E3. 
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MU > 0 
Edge 2 


Node-Line Offset for Node-Line Offset for 

Plate Strip No. 7 Plate Strip No. 9 


Figure E3. Unit-cell cross-sectional representation of a hat-stiffened panel, 
including the geometric details of the node-line offsets. 
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File ’pasco.in* for a hat-stiffened panel. 


HAT - UNIT-CELL REPRESENTATION OF A HAT-STIFFENED PANEL 

&CONDAT 

&END 

&PANEL 

B=2.,.75, 2.267855,. 75, 2.,.75, 1.5,1., 1.5, .75, 

EL=15„ 

CONV1=50000„ 

FREQ=0„ 

IP=0, 

MAXJJJ=0, 

MAT(1,1)=1,1,1, 

T(l)=.005,.005,.005,.07,.0634,.0296, 

THET(I)=0, 45,90, 

K WALL( 1 , 1 )=2,-2,2,-2, 1,3, 

KWALL( 1 ,2)=2,-2, 1 ,3,2, -2, 1 ,3, 

KWALL(1,3)=2,-2,1,3,2,-2,1,3, 

KWALL(1 ,4)=2,-2, 1 ,3, 1 , 1 ,2,-2, 1 ,3, 

IWALL(1)=1,1,1,1,1,2,3,4,3,2, 

FSTIFF=10„ 

HCARD=6,-1 1,-6, 0,-4 ,0,-4, 

6,- 12,-7, -5, -6,0,0, 

6,-1 3,-9, 0,0,5,-6, 

6,-14,-10,0,-4,0,-4, 

4,-15,12,-65,0, 

4,-16,13,65,0, 

4,17,15,8,16, 

14,18,1 ,-990,9000,2,- 1 1 ,3,- 1 7,4,- 1 4,5,0,-990,9000, 

ICARD=5, 1,2, 1,-990.9000, 

5,2,3,2,3,11, 

3,3,6, 3, 

3.3.4.15, 

3.4.5. 8, 

3.5.6. 16, 

5,6,7,4,7,14, 

3,7,8,5, 

3.8, -990,9000, 

ICREP=1, 

NOBAY=l, 

NX(1)=1„ 

NY(1)=0., 

&END 

&MATER 

E1(1)=18.5E6, E2(1)=1.64E6, E12(1)=.87E6, ANU1(1)=.30, 

RHO( 1)= 1 .477E-4, ALFA 1 ( 1 )=0.25E-6, ALFA2(l)=16.2E-6, 
ALLOW(1,1)=2,.01 1,-0.011,0.01,-0.013,0.0155. 

&END 

— End of file ’pasco.in’ — 
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File ’nlpan.inl’ for a hat-stiffened panel. 

HAT - UNIT-CELL REPRESENTATION OF A HAT-STIFFENED PANEL 


0000 
0000 
’STOLL’ 1 0 
’Y’ 1 8 

2 ’Y’ 

1 2 

1 5 
1 

1 3 
10 8 
4 

60-40-4 
7-5-600 
9005-6 
10 0-4 0-4 
0 . 1 1 0 . 

0 . 1 1 0 . 

0 . 1 1 0 . 

0 . 1 1 0 . 

0 . 1 1 0 . 

0.0 0 0 . 

-65. 0 0 0. 

0.0 0 0 . 

65. 0 0 0. 

0.0 0 0 . 

1 2 
23 

3 6 
67 
7 8 
34 

4 5 

5 6 
67 
1 8 

2 3 2 1 
2 3 2 1 
00 
00 


! ICOORP JLNPRT.IEFPRT.IQUPRT 
! INGPRT,ICOPRT,IPRINT t IPBUG 
IUSER, IHALT, MSRCH 
ISYMSTR, NODSYM(l), NODSYM(2) 

INMUSE, FORCE 
!NSOL(I), I=1,NMUSE 
!MHIN(I), I=1,NMUSE 
!ISOL(l,J), J=l,NSOL(I), 1=1 
! 1=2 
INPLATE NNODE 
INPOFFS 

!IP,IECVY(IP,1)JECVZ(IP,1)JECVY(IP^),IECVZ(IP^) IPOFFS 
I IPOFFS =2 

I IPOFFS=3 

I IPOFFS =4 

IMU(IP), SKVEC(IP), PRVEC(IP) TUA(IP), IP=1 

IP=2 
IP=3 
IP=4 
IP=5 
IP=6 
IP=7 
IP=8 
IP=9 
IP=10 

!NOD(IP,IE), IP=1 

EP=2 
IP=3 
IP=4 
IP=5 
IP=6 
IP=7 
IP=8 
IP=9 

BNODE(IB), IB=1,2 
BCVEC(IC,IB), IC=1 to 4, IB=1 
IB=2 

NCU, NCW 
NSPU, NSPW 


**** 


**** DESCRETIZATION AND REFERENCE VALUES 
12 ININTN 

1 0 3 0 IIPHREF, IPWDET, IPDFL, IYDFL 

**** LOADING AND MODELLING OPTIONS **** 

’D’ ’N’ ICONTRL, COMPAT 

1 1. 0. 0. 0. I ILSET .NXGLJVY GL,EPXL,EPYLG 

0.0 IPRUNIT 

’Y’ ILOCGLO 

2 IINDPLT 

**** LOAD, WIDTH, WALL INDEXING **** 
1232145654 !ILBWAL(IP),IP=1,NPLATE 
4 INKWALL 

-- End of file ’nlpan.inl’ — 
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E2. Stoll, F., "Geometrically Nonlinear Analysis of Stiffened Composite Panels with Various End-Support 
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Structural Dynamics, and Materials Conference, La Jolla, CA, April 19-21, 1993. 
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Type of 
Load Control 


Load Axis 1 


Load Axis 2 


Option 1 Option 2 


( BCVEC(2,IB)=1, ( BCVEC(2,IB)=2 or 3, 
IB=1,2 ) IB=1 or IB=2 ) 


Displacement 

Control 

(CONTRL=’D') 

Aw = XAu l 

Av = \Av l 

o 

II 

£ 

Load 

Control 

(CONTRL='L') 


o 

II 

\< 

NyG = MyG L 


Table I. Options for Control of the In-Plane Loading. IB = boundary index number (1 or 2). 
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Table 2. Input Parameters Corresponding to Several Different Cases for In-Plane Loading. 
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Component 

(IC) 

Option 1 

(BCVEC(IC,IB)=1) 

Option 2 

(B CVEC(I C,IB)=2) 

Option 3 

(BCVEC(IC,IB)=3) 

1: U\F; 

u n lxx =o 

f;=o 

- 

HH 

V 1 

II 

i 

II 

v,:=o 

3: W\ F; 

W n = 0 

f;=o 

- 

4: V F", M " 

o 

II 

M n = 0 

- 


Table 3. Options Available for Specifying Conditions Along the Boundary Node-Lines. IB = boundary index 
number (1 or 2). 


Case 

Simplification ] 

CONTRL = D’ 

/"S 

ri 

II 

*»» 

o 

III 


£ 

III 

o 

O 

II 

£ 

F 

in 

o 

-or- 


BCVEC(2,IB)=1, 

(IB=1,2) 


BCVEC(2,IB) * 3 

J30 

III 

o 

/— \ 

fN # 

II 

JO 

W' 

(IB=1 or 2) 

{u fi } so 


Table 4. Simplifications to the General Displacement Form for Special Cases. 
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Local 

Postbuckling 0 


Global 

Postbuckling, 
s.s. ends 


Global 

Postbuckling, 
clamped ends 


Local/Global 
Mode Interaction, 
s.s. ends 


Local/Global 
Mode Interaction, 
clamped ends 


Symmetrical Structure 


1. (m,, 1) 

2. (m,,i 2 ) 

(3m,, »',) (3 m,,/ 2 ) 

3. (m,,i 3 ) 

(5 m,,/,) 


1 . ( 1 , 1 ) 

2. (3^) (3a) 


1 . ( 1 , 1 ) 

(3a) (3a) 

2. (5a) (5 a) 


Al. (1,1) 

(rru,U\) (i m,,ui ) 

2. (jn&2,U\) (w2,m 2 ) 
(nu-2,Ui) (mr2,u 2 ) 

Bl. (1,1) 

(wwi) (w,a) 

2. (/nrh2,j,) (mtf-2A) 

(/Ht2a) (wir2,A) 


Al. (1.1) 

(3a) (3,50 

(m,,ui) (m,,« 2 ) 

(/th-2,«i) (mri-2,u 2 ) 

(jru-2,U\) ( mr2,u 2 ) 


Unsymmetrical Structure 


1. (mu.l) (m,,2)* ... 

2. (3m,,l) (3m,,2) ... 

3. (5m,. 1 ) (5m,, 2) ... 


1. (1,1) (1,2)* 

2. (3,1) (3,2) ... 


1. G.l) (1,2)* 
(3,1) (3,2) ... 

2. (5,1) (5,2) ... 


1 . ( 1 , 1 ) ( 1 . 2 )* 

(m<,l) (m,,2) ... 

2. (mn-2,1) (mn-2,2) ... 

(to- 2,1) (mr2,2) ... 


1. (1,1) 0,2)* 

(3,1) (3,2) ... 

(m,,l) (m,, 2) ... 
(mn-2,1) (/*+ 2,2) ... 
(mr2,l) (mr2,2) ... 


Bl. (1,1) 

(3a) (3a) 

(m,A) (m,A) 

(, m & 2a) (mr(-2,5 2 ) 
(mr2A) (mr2 a) 


m, - Critical halfwave number for local buckling. 

(m,«) - Unsymmetric mode. 

(m,s) - Symmetric mode. 

1. - Essential modes. 

2., 3. - Supplemental modes for improved accuracy. 

A, B - Dual mode sets which should each be used independently in separate analyses. 

* Modes (m.i,). (m,/ 2 ), ... match mode (m,,l) w.r.t. transverse symmetry or lack thereof. 
‘ Include mode (m,2) if its eigenvalue is close to that for (m,l). 
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Figure 1. Generalized load types modelled by NLPAN. 
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Figure 2. Labeling conventions for a representative plate strip. 






a) Symbol definitions. 



Figure 4. Relative orientation of the side -edge of a plate strip and its associated node line. 
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3.05 cm 


50 plies 


22 plies 


3.18 cm 


12 plies 


4 plies 


\ 


i i — 16 plies 

i [t 45/0 2 / +45/90 2 ] s 


— 4.32 cm 

b) Stiffener details. 



c) Profile of the theoretical buckling mode for a unit cell representation with 17 plate strips 
and 12 node lines (buckling mode has 5 longitudinal halfwaves). 

Figure 5. 1-stiffened graphite/epoxy panel configuration. 
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Extension of ..•** o'' 

Prebuckling Path 



0 12 3 4 

Au/ Auer 


a) Normalized end load versus normalized end shortening. 



Figure 6. Analytical and experimental results for an I-stiffened panel. 
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Gage C 


age A 


''-Gage D M3age B 

a) Strain gage placement, mid-length, center bay. 




c) Strains at a stiffener flange. 


Figure 7. Longitudinal surface strains on an I-stiffened panel. 
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a) Unit-stiffener-cell analysis model. 


p limit /Pcr E 



Normalized Bowing Imperfection, w°/H 

b) Limit loads from experiments, and from analysis with and without 
orthogonality imposed between {ujj} and {uj }. 

Figure 8. Results for an imperfection-sensitive thin-blade-stiffened isotropic panel 
loaded in uniaxial compression. 
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Boundary conditions common to ail analyses: 

x=-10, 10: u=0, N X y=0, w=0 
y=-5, 5: u=0, v=0, w=0, w, y =0 



a) Test section and boundary conditions 



b) Stiffener and laminate details. 


Figure 9. Pressure-loaded panel configuration. 
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w/h w/h 


8 


6 

i 

2 

0 


-0.5 -0.25 0 0.25 0.5 

y/B 

a) Displacements across the width at the mid-length. 



b) Displacements along the length under the stiffener (y=0). 



x/L 

c) Displacements along the length on the skin (y/B=. 25) 
Figure 10. Skin displacements with one atmosphere pressure (14.7 psi). 
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Lamina properties: 

El = 22.5 Msi 
E2 = 1.17 Msi 
G12 = 0.66 Msi 
v 12 = 0.22 

a 1 = -.04* 10"® per degree F 
a 2 = 16.7*10'® per degree F 
h p|y = 0.005 in. 

Laminate configuration: [+45/-45/0/0]s 

Boundary conditions: x=0, L: u=0, w=0, M x =0, N xy =0 

y=0, B: v=0, w=0, M y =0. N xy =0 

Buckling temperature: Tcrit = 69.4 deg. F (Ref. [28]) 

Tcrit = 71 .4 deg. F (NLPAN) 



a) Plate configuration and properties. 



Figure 11. Laminated composite plate subjected to thermal loading. 
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